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

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

      // CM TopTagger (unpublished)
      CMTopTagger cmtagger(cs, jet, zmin);
      if (cmtagger.maybe_top(top_range, W_range)) {
	cout << "      CM tagged       " 
	     << "mt = " << cmtagger.top_subjet().m()
	     << " mW = " << cmtagger.W_subjet().m() << endl;
      }

    }

    
  }
}