#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/SISConePlugin.hh"
#include "UWEvent.hh"
#include "print_jets.hh"
#include "SimpleHist.hh"

using namespace std;
using namespace fastjet;


//
// Run this program with
//    gunzip -c ../event-files/pythia-dijets.UW.gz | ./02-jet-definitions
//

int main() {

  // max number of events read in (even if input file contains more)
  int maxiev = 10000;
  
  // filename containing the events
//  string datafile =  "./data/";
//  datafile += "pythia-dijets.UW";
  
  // define a vector of many jet definitions
  vector<JetDefinition> jet_defs;
  // antikt, R=0.6
  jet_defs.push_back(JetDefinition(antikt_algorithm,0.6));
  // antikt, R=0.4
  jet_defs.push_back(JetDefinition(antikt_algorithm,0.4));
  // kt, R=0.4
  jet_defs.push_back(JetDefinition(kt_algorithm,0.4));
  // SISCONE, R=0.4, f=0.75
  double overlap_threshold = 0.75;
  SISConePlugin siscone(0.4,overlap_threshold);
  jet_defs.push_back(JetDefinition(&siscone));

  // initialize histograms
  vector<SimpleHist> inclusive_jets_hist(jet_defs.size(),SimpleHist(0.0,100.0,1.0));  
  
  // read file containing events
//  ifstream input_file;
//  input_file.open(datafile.c_str());
  int iev=0;
  vector<PseudoJet> input_particles;
  while (readUWEvent(cin, input_particles) && ++iev <= maxiev) {
//  while (readUWEvent(input_file, input_particles) && ++iev <= maxiev) {
    cerr << "Event " << iev << ", size = " << input_particles.size() << endl; 

    // run over various jet definitions
    for (unsigned int j=0; j < jet_defs.size(); j++) {
       // cluster the input particles
       ClusterSequence cs(input_particles, jet_defs[j]);
  
       // extract the jets
       vector<PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());

       // add ALL jets to histogram
       for (unsigned int i=0; i < jets.size(); i++) {
          inclusive_jets_hist[j].add_entry(jets[i].perp());
       }
    
       // print out jets in first event
       if (iev < 1 ) { print_jets(jets); } // in print_jets.hh
   }
  }
  // open output file for histograms
  string outfile = "02-jet-definitions.hist";
  ofstream output_file;
  output_file.open(outfile.c_str());
  // run over various jet definitions
  for (unsigned int j=0; j < jet_defs.size(); j++) {
      // output histograms
      output_file << "# j = " << j << endl;
      output_file << "# " << jet_defs[j].description() << endl;
      output(inclusive_jets_hist[j],&output_file); // in SimpleHist.hh
      output_file << endl << endl;
  }
  output_file.close();
  
  
}