#ifndef  D0RunIIconeJets_CONESPLITMERGE
#define  D0RunIIconeJets_CONESPLITMERGE
// ---------------------------------------------------------------------------
// ConeSplitMerge.hpp
//
// Created: 28-JUL-2000 Francois Touze
//
// Purpose: Implements the pT ordered jet split/merge algorithm for the 
//   Improved Legacy Cone Algorithm split/merge algo.
//
// Modified:
//   31-JUL-2000  Laurent Duflot
//     + introduce support for additional informations (ConeJetInfo)
//    1-AUG-2000  Laurent Duflot
//     + jet merge criteria was wrong, now calculate shared_ET and compare to 
//       neighbour jet pT.
//     + split was incomplete: shared items not really removed from jets.
//    4-Aug-2000  Laurent Duflot
//     + use item methods y() and phi() rather than p4vec() and then P2y and 
//       P2phi
//    7-Aug-2000  Laurent Duflot 
//     + force the list to be organized by decreasing ET and, for identical ET,
//       by decreasing seedET. Identical protojets can be found by eg nearby
//       seeds. The seedET ordering is such that for identical jets, the one
//       with the highest seedET is kept, which is what we want for efficiency
//       calculations.
//     + avoid unecessary copies of lists by using reference when possible
//    9-Aug-2000  Laurent Duflot
//     + save initial jet ET before split/merge
//    1-May-2007 Lars Sonnenschein
//    extracted from D0 software framework and modified to remove subsequent dependencies
//
//
// This file is distributed with FastJet under the terms of the GNU
// General Public License (v2). Permission to do so has been granted
// by Lars Sonnenschein and the D0 collaboration (see COPYING for
// details)
//
// History of changes in FastJet compared tothe original version of
// ConeSplitMerge.hpp
//
// 2011-12-13  Gregory Soyez  <soyez@fastjet.fr>
// 
//        * added license information
//
// 2009-01-17  Gregory Soyez  <soyez@fastjet.fr>
//
//        * put the code in the fastjet::d0 namespace
//
// 2007-12-14  Gavin Salam  <salam@lpthe.jussieu.fr>
// 
//        * replaced make_pair by std::make_pair
//
// ---------------------------------------------------------------------------


#include <iostream>
#include <map>
#include <utility>
#include <vector>
#include "ProtoJet.hpp"

//using namespace D0RunIIconeJets_CONEJETINFO;

#include <fastjet/internal/base.hh>

FASTJET_BEGIN_NAMESPACE

namespace d0{

//
// this class is used to order ProtoJets by decreasing ET and seed ET
template <class Item>
class ProtoJet_ET_seedET_order
{
public:
  bool operator()(const ProtoJet<Item> & first, const ProtoJet<Item> & second) const
  {
    if ( first.pT() > second.pT() ) return true;
    else
      if ( first.pT() < second.pT() ) return false;
      else return ( first.info().seedET() > second.info().seedET() );
  }
};


template <class Item>
class ConeSplitMerge {

public :

  typedef std::multimap<ProtoJet<Item>,float,ProtoJet_ET_seedET_order<Item> > PJMMAP;
  
  ConeSplitMerge();
  ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector);
  ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist);
  ~ConeSplitMerge() {;}
  void split_merge(std::vector<ProtoJet<Item> >& ecv,float s, float pT_min_leading_protojet, float pT_min_second_protojet, int MergeMax, float pT_min_noMergeMax);

private :
  PJMMAP _members;
};
///////////////////////////////////////////////////////////////////////////////
template<class Item>
ConeSplitMerge<Item>::ConeSplitMerge():_members() {;}

template<class Item>
ConeSplitMerge<Item>::ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector) 
{
  // sort proto_jets in Et descending order
  typename std::vector<ProtoJet<Item> >::const_iterator jt;
  for(jt = jvector.begin(); jt != jvector.end(); ++jt) 
  {
    // this is supposed to be a stable cone, declare so
    ProtoJet<Item> jet(*jt);
    jet.NowStable();
    _members.insert(std::make_pair(jet,jet.pT()));
  }
}

template<class Item>
ConeSplitMerge<Item>::ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist) 
{
  //_max_nb_items =-1;
  // sort proto_jets in Et descending order
  typename std::list<ProtoJet<Item> >::const_iterator jt;
  for(jt = jlist.begin(); jt != jlist.end(); ++jt) 
  {
    // this is supposed to be a stable cone, declare so
    ProtoJet<Item> jet(*jt);
    jet.NowStable();
    _members.insert(std::make_pair(jet,jet.pT()));
  }
}

template<class Item>
void ConeSplitMerge<Item>::split_merge(std::vector<ProtoJet<Item> >& jcv,
				       float shared_ET_fraction,
				       float pT_min_leading_protojet, float pT_min_second_protojet,
				       int MergeMax, float pT_min_noMergeMax) 
{
  while(!_members.empty()) 
  {
    /*
    {
      std::cout << std::endl;
      std::cout << " ---------------  list of protojets ------------------ " <<std::endl;
      for ( PJMMAP::iterator it = _members.begin();
	    it != _members.end(); ++it)
      {
	std::cout << " pT y phi " << (*it).pT() << " " << (*it).y() << " " << (*it).phi() << " " << (*it).info().seedET() <<  " " << (*it).info().nbMerge() << " " << (*it).info().nbSplit() << std::endl;
      }
      std::cout << " ----------------------------------------------------- " <<std::endl;
    }
    */


    // select highest Et jet
    typename PJMMAP::iterator itmax= _members.begin();
    ProtoJet<Item> imax((*itmax).first);
    const std::list<const Item*>& ilist(imax.LItems());

    _members.erase(itmax);
 
    // does jet share items?
    bool share= false;
    float shared_ET = 0.;
    typename PJMMAP::iterator jtmax;
    typename PJMMAP::iterator jt;
    for(jt = _members.begin(); jt != _members.end(); ++jt) 
    {
      const std::list<const Item*>& jlist((*jt).first.LItems());
      typename std::list<const Item*>::const_iterator tk;
      int count;
      for(tk = ilist.begin(), count = 0; tk != ilist.end(); 
	  ++tk, ++count) 
      {
	typename std::list<const Item*>::const_iterator where= 
	  find(jlist.begin(),jlist.end(),*tk);   
	if(where != jlist.end()) 
	{
	  share= true;
	  shared_ET += (*tk)->pT();
	}
      }
      if(share) 
      {
	jtmax = jt;
	break;
      }
    }
    
    if(!share) {
      // add jet to the final jet list
      jcv.push_back(imax);
      //std::cout << " final jet " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 
    }
    else 
    {
      // find highest Et neighbor
      ProtoJet<Item> jmax((*jtmax).first);

      // drop neighbor
      _members.erase(jtmax);


      //std::cout << " split or merge ? " << imax.pT() << " " << shared_ET << " " << jmax.pT() << " " << s << " " << (jmax.pT())*s << std::endl;

      // merge
      if( shared_ET > (jmax.pT())*shared_ET_fraction 
	  && (imax.pT() > pT_min_leading_protojet || jmax.pT() > pT_min_second_protojet)
	  && (imax.info().nbMerge() < MergeMax || imax.pT() > pT_min_noMergeMax))
	{
	  // add neighbor's items to imax
	  const std::list<const Item*>& jlist(jmax.LItems());
	  typename std::list<const Item*>::const_iterator tk;
	  typename std::list<const Item*>::const_iterator iend= ilist.end();
	  bool same = true; // is jmax just the same as imax ? 
	  for(tk = jlist.begin(); tk != jlist.end(); ++tk) 
	    {
	      typename std::list<const Item*>::const_iterator where= 
		find(ilist.begin(),iend,*tk);   
	      if(where == iend) 
		{
		  imax.addItem(*tk);
		  same = false;
		}
	    }
	  if ( ! same ) 
	    {
	      // recalculate
	      //float old_pT = imax.pT();
	      
	      imax.updateJet();
	      imax.merged();
	      //std::cout << " jet merge :: " << old_pT << " " << jmax.pT() << " " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 
	    }
	}
      
      //split and assign removed shared cells from lowest pT protojet
      else if(shared_ET > (jmax.pT())*shared_ET_fraction)
      {
	// here we need to pull the lists, because there are items to remove                                                                           
	std::list<const Item*> ilist(imax.LItems());
	std::list<const Item*> jlist(jmax.LItems());

        typename std::list<const Item*>::iterator tk;
        for(tk = jlist.begin(); tk != jlist.end(); )
	  {
	    typename std::list<const Item*>::iterator where=
	      find(ilist.begin(),ilist.end(),*tk);
	    if(where != ilist.end()) {
	      tk = jlist.erase(tk);
	    }
	    else ++tk;
	  }
	
        jmax.erase();

        for ( typename std::list<const Item*>::const_iterator it = jlist.begin();
              it != jlist.end(); ++it) jmax.addItem(*it);

        // recalculated jet quantities 
        jmax.updateJet();
        jmax.splitted();
        //std::cout << " jet split 1 :: " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl;                         
        _members.insert(std::make_pair(jmax,jmax.pT()));
      }

      // split and assign shared cells to nearest protojet
      else 
      {
	// here we need to pull the lists, because there are items to remove
	std::list<const Item*> ilist(imax.LItems());
	std::list<const Item*> jlist(jmax.LItems());
	
	typename std::list<const Item*>::iterator tk;
	for(tk = jlist.begin(); tk != jlist.end(); ) 
	{
	  typename std::list<const Item*>::iterator where= 
	    find(ilist.begin(),ilist.end(),*tk);
	  if(where != ilist.end()) {
	    float yk   = (*tk)->y();
	    float phik = (*tk)->phi();
	    float di= RD2(imax.y(),imax.phi(),yk,phik);
	    float dj= RD2(jmax.y(),jmax.phi(),yk,phik);
	    if(dj > di) 
	    {
	      tk = jlist.erase(tk);
	      //std::cout << " shared item " << (*tk)->pT() << " removed from neighbour jet " << std::endl;
	    }
	    else
	    {
	      ilist.erase(where);
	      ++tk;
	      //std::cout << " shared item " << (*tk)->pT() << " removed from leading jet " << std::endl;
	    }
	  }
	  else ++tk;
	}
	// recalculate jets imax and jmax
	
	// first erase all items
	imax.erase();
	// put items that are left
	for ( typename std::list<const Item*>::const_iterator it = ilist.begin();
	      it != ilist.end(); ++it) imax.addItem(*it);
	// recalculated jet quantities
	imax.updateJet();
	imax.splitted();
	//std::cout << " jet split 2 :: " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 


	// first erase all items
	jmax.erase();
	// put items that are left
	for ( typename std::list<const Item*>::const_iterator it = jlist.begin();
	      it != jlist.end(); ++it) jmax.addItem(*it);
	// recalculated jet quantities
	jmax.updateJet();
	jmax.splitted();
	//std::cout << " jet split " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl; 

	_members.insert(std::make_pair(jmax,jmax.pT()));
      }
      _members.insert(std::make_pair(imax,imax.pT()));
    }
  } // while
}
///////////////////////////////////////////////////////////////////////////////

}  // namespace d0

FASTJET_END_NAMESPACE

#endif