#include "fastjet/ClusterSequence.hh"
#include "UWEvent.hh"
#include "JHTopTagger.hh"
#include <iostream>

using namespace std;
using namespace fastjet;

int main() {
  double R = 0.6;
  JetDefinition jet_def(cambridge_algorithm, R);

  vector<PseudoJet> event;
  int iev = 0;

  while (readUWEvent(cin, event) && ++iev <= 5) {
    cout << "Event " << iev << ", size = " << event.size() << endl;

    // run the jet clustering
    ClusterSequence cs(event, jet_def);
    double ptmin = 500.0;
    vector<PseudoJet> jets = sorted_by_pt(cs.inclusive_jets(ptmin));

    // loop over the jets
    for (unsigned i = 0; i < jets.size(); i++) {
      const PseudoJet & jet = jets[i];
      cout << "  jet " << i << ":" 
	   << "pt = " << jet.perp() << ", "
	   << "phi = " << jet.phi() << ", "
	   << "m = "  << jet.m() << endl;
      
      // JH TopTagger parameters 
      double zmin = 0.05; // momentum fraction cut
      double drmin = 0.19; // distance cut
      Range  top_range(145, 205), W_range(65, 95);
      double max_cos_theta_h = 0.7;
      JHTopTagger jhtagger(cs, jet, zmin, drmin);
      if (jhtagger.maybe_top(top_range, W_range, max_cos_theta_h)) {
	cout << "      JH tagged       " 
	     << "mt = " << jhtagger.top_subjet().m()
	     << " mW = " << jhtagger.W_subjet().m() << endl;
      }

      // run Ysplitter
      vector<PseudoJet> constituents = cs.constituents(jet);
      JetDefinition ktdef(kt_algorithm, 1.0);
      ClusterSequence cskt(constituents, ktdef);
      PseudoJet ktjet = sorted_by_pt(cskt.inclusive_jets())[0];
      double dmerge12 = cskt.exclusive_subdmerge(ktjet, 1); 
      double dmerge23 = cskt.exclusive_subdmerge(ktjet, 2); 
      cout << "           sqrt(dmerge) " << sqrt(dmerge12) << " " << sqrt(dmerge23) << endl;
      
    }

    
  }
}