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

using namespace std;
using namespace fastjet;


//
// Run this program with
//    gunzip -c ../event-files/pythia-dijets-100-110.UW.gz | ./03-areas
//

int main() {

  // select the kind of jets to look at
  //
  // NB. It is not efficient to do these one at a time, because we recluster the
  // event each time. Ideally, we'd extract and histogram all quantities of interest
  // in a single run.
  double rapmax = 2.;
  Selector sel_rapmax = SelectorAbsRapMax(rapmax);
  // in each event select only 2 hardest jets, followed by selection of jets in |y| < 2
  Selector full_selector = sel_rapmax * SelectorNHardest(2);
  string outfile = "03-areas-hard.hist";
  // select all BUT the 4 hardest jets in each event, excluding the pure-ghost jets
//  Selector full_selector = ((!SelectorNHardest(4))&&(!SelectorIsPureGhost())) * sel_rapmax;
//  string outfile = "03-areas-soft.hist";
  // select only pure-ghost jets
//  string outfile = "03-areas-ghost.hist";
//  Selector full_selector = SelectorIsPureGhost() * sel_rapmax;

  // max number of events read in (even if input file contains more)
  int maxiev = 1000;
  
  // 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;                 //  VERY slow with areas
//  SISConePlugin siscone(0.4,overlap_threshold);    //  Run this at your own peril
//  jet_defs.push_back(JetDefinition(&siscone));     //

  // set up area definition
  AreaDefinition area_def(active_area_explicit_ghosts,GhostedAreaSpec(rapmax + 1.));

  // initialize histograms
  vector<SimpleHist> inclusive_jets_hist(jet_defs.size(),SimpleHist(0.0,3.,0.02));  

  // initialize AverageAndError
  vector<AverageAndError> avg_area(jet_defs.size());
  vector<AverageAndError> avg_pt(jet_defs.size());
  
  // 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
       ClusterSequenceArea cs(input_particles, jet_defs[j], area_def);
  
       // extract the jets
       vector<PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());

       // select the jets
       jets = full_selector(jets);

       // add all selected jets to histogram, add entries to averages
       for (unsigned int i=0; i < jets.size(); i++) {
           double area_over_piR2 = jets[i].area()/pi/pow(jet_defs[j].R(),2);
           inclusive_jets_hist[j].add_entry(area_over_piR2);
           avg_area[j].add(area_over_piR2);
           avg_pt[j].add(jets[i].perp());
       }
    
       // print out jets in first event
       //if (iev < 1 ) { print_jets(jets); } // in print_jets.hh
   }
  }
  // open output file for histograms
  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_file << "# " << area_def.description() << endl;
      output_file << "# jets selected using: " << full_selector.description() << endl;
      output_file << "# average area = " << avg_area[j].average() << endl;
      output_file << "# average pt   = " << avg_pt[j].average() << endl;
      output(inclusive_jets_hist[j],&output_file); // in SimpleHist.hh
      output_file << endl << endl;
  }
  output_file.close();
  
  
}