// $Id$
//
// Copyright (c) 2013,
// Andrew Larkoski, Simone Marzani, Gregory Soyez, and Jesse Thaler
//
// Run this exanmple with 
//
//   example < ../data/single-event.dat
//
//----------------------------------------------------------------------
// This file is part of FastJet contrib.
//
// It is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at
// your option) any later version.
//
// It is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
// License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this code. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------

#include <iostream>
#include <sstream>

#include "fastjet/ClusterSequence.hh"
#include "fastjet/Selector.hh"

#include "SoftDrop.hh" // In external code, this should be fastjet/contrib/SoftDrop.hh

using namespace std;
using namespace fastjet;

// forward declaration to make things clearer
void read_event(vector<PseudoJet> &event);

//----------------------------------------------------------------------
int main(){

  //----------------------------------------------------------
  // read in input particles
  vector<PseudoJet> event;
  read_event(event);
  cout << "# read an event with " << event.size() << " particles" << endl;

  //----------------------------------------------------------
  // illustrate how this SoftDrop contrib works
  //
  // The plan goes as follows:
  //  - cluster the event and get the 2 hardest jets
  //  - for each of these 2 jets, apply the soft drop and show a few
  //    properties of the resulting jet

  // Cluster the event and get the initial 2 hardest jets
  JetDefinition jet_def(antikt_algorithm, 0.6);
  ClusterSequence cs(event, jet_def);
  vector<PseudoJet> jets = SelectorNHardest(2)(cs.inclusive_jets());

  // declare a SoftDropTagger to play with
  double beta = 1.0;
  double zcut = 0.1;
  double mu   = 1.0;
  contrib::SoftDropTagger soft_drop(beta, zcut, mu);

  // apply the soft drop to the jets
  vector<PseudoJet> sd_jets = soft_drop(jets);

  // show some info for all these jets
  for (unsigned int i=0; i<jets.size(); i++){
    // create handful shortcuts for te original and soft-dropped jets
    const PseudoJet &jet    = jets[i];
    const PseudoJet &sd_jet = sd_jets[i];

    // print the jet kinematic properties:
    cout << "Original jet: pt=" << jet.pt() 
         << ", y=" << jet.rap() 
         << ", phi=" << jet.phi()
         << ", mass=" << jet.m() << endl;

    // if soft-drop failed, just skip all the information
    if (sd_jet==0){
      cout << "Soft-dropped jet: zero PseudoJet" << endl;
      continue;
    }

    cout << "Soft-dropped jet: pt=" << sd_jet.pt() 
         << ", y=" << sd_jet.rap() 
         << ", phi=" << sd_jet.phi() 
         << ", mass=" << sd_jet.m() << endl;

    // show the last splitting information
    cout << "soft-drop value = " << sd_jet.structure_of<contrib::SoftDropTagger>().soft_drop()
         << endl
         << "mass-drop value = " << sd_jet.structure_of<contrib::SoftDropTagger>().mu()
         << endl;

    // get the 2 prongs in the tagged jet
    cout << "Distance between the 2 prongs=" << sd_jet.structure_of<contrib::SoftDropTagger>().Rg() << endl;

    // show the energy loss
    cout << "Energy loss = " << 1-sd_jet.pt()/jet.pt() << endl;

    // show the z_drop_max
    // (Commented out because still in beta testing)
    cout << "# z_drop_max = " << sd_jet.structure_of<contrib::SoftDropTagger>().z_drop_max() << endl;

    

    cout << endl;
  }

  return 0;
}

// read in input particles
void read_event(vector<PseudoJet> &event){  
  string line;
  while (getline(cin, line)) {
    istringstream linestream(line);
    // take substrings to avoid problems when there are extra "pollution"
    // characters (e.g. line-feed).
    if (line.substr(0,4) == "#END") {return;}
    if (line.substr(0,1) == "#") {continue;}
    double px,py,pz,E;
    linestream >> px >> py >> pz >> E;
    PseudoJet particle(px,py,pz,E);

    // push event onto back of full_event vector
    event.push_back(particle);
  }
}