/*
 *   ---- Delphes ----
 *     A Fast Simulator for general purpose LHC detector
 *       S. Ovyn ~~~~ severine.ovyn@uclouvain.be
 *
 *         Center for Particle Physics and Phenomenology (CP3)
 *           Universite Catholique de Louvain (UCL)
 *             Louvain-la-Neuve, Belgium
 *             */

// \brief Trigger class, and some generic definitions

#include "JetsUtil.h"

using namespace std;

JetsUtil::JetsUtil() : plugins(NULL) {
   DET = new RESOLution();
   init();
}


JetsUtil::JetsUtil(const string & DetDatacard) : plugins(NULL) {
   DET = new RESOLution();
   DET->ReadDataCard(DetDatacard);
   init();
}


JetsUtil::JetsUtil(const RESOLution * DetDatacard) : plugins(NULL) {
   DET = new RESOLution(*DetDatacard);
   init();
}


JetsUtil::JetsUtil(const JetsUtil& jet){
   DET = new RESOLution(*(jet.DET));
}

JetsUtil& JetsUtil::operator=(const JetsUtil& jet) {
   if(this==&jet) return *this;
   DET = new RESOLution(*(jet.DET));	
   return *this;
}

void JetsUtil::init() {
  if(plugins) delete plugins;
 
  switch(DET->JET_jetalgo) {
  default:
  case 1: {
    // set up a CDF midpoint jet definition
#ifdef ENABLE_PLUGIN_CDFCONES
    plugins = new fastjet::CDFJetCluPlugin(DET->JET_seed,DET->JET_coneradius,DET->JET_C_adjacencycut,DET->JET_C_maxiterations,DET->JET_C_iratch,DET->JET_overlap);
    jet_def = fastjet::JetDefinition(plugins);
#else
    plugins = NULL;
#endif
  }
    break;
    
  case 2: {
    // set up a CDF midpoint jet definition
#ifdef ENABLE_PLUGIN_CDFCONES
    plugins = new fastjet::CDFMidPointPlugin (DET->JET_seed,DET->JET_coneradius,DET->JET_M_coneareafraction,DET->JET_M_maxpairsize,DET->JET_M_maxiterations,DET->JET_overlap);
    jet_def = fastjet::JetDefinition(plugins);
#else
    plugins = NULL;
#endif
  }
    break;
    
  case 3: {
    // set up a siscone jet definition
#ifdef ENABLE_PLUGIN_SISCONE
    plugins = new fastjet::SISConePlugin (DET->JET_coneradius,DET->JET_overlap,DET->JET_S_npass, DET->JET_S_protojet_ptmin);
    jet_def = fastjet::JetDefinition(plugins);
#else
    plugins = NULL;
#endif
  }
    break;
    
  case 4: {

    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->JET_coneradius);
  }
    break;
    
  case 5: {
    jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,DET->JET_coneradius);
  }
    break;
    
  case 6: {
    jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,DET->JET_coneradius);
  }
    break;
  }
  
}

vector<fastjet::PseudoJet> JetsUtil::RunJets(const vector<fastjet::PseudoJet>&  input_particles)
{
  inclusive_jets.clear();
  sorted_jets.clear();
  // run the jet clustering with the above jet definition
  if(input_particles.size()!=0)
    {
      fastjet::ClusterSequence clust_seq(input_particles, jet_def);
      // extract the inclusive jets with pt > 5 GeV
      double ptmin = 5.0;
      inclusive_jets = clust_seq.inclusive_jets(ptmin);
      // sort jets into increasing pt
      sorted_jets = sorted_by_pt(inclusive_jets);
    }
  
  return sorted_jets;                                                    
}

void JetsUtil::RunJetBtagging(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchJet,const vector<fastjet::PseudoJet> & sorted_jets,const TSimpleArray<TRootGenParticle>& NFCentralQ)
{
  TRootJet *elementJet;
  TLorentzVector JET;
  for (unsigned int i = 0; i < sorted_jets.size(); i++) {
    JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
    if(JET.Pt() > DET->PTCUT_jet)
      {
	elementJet = (TRootJet*) branchJet->NewEntry();
	elementJet->Set(JET);
	// b-jets
	bool btag=false;
	if((fabs(JET.Eta()) < DET->CEN_max_tracker && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
	elementJet->Btag = btag;
      }
  } // for itJet : loop on all jets
  
}

void JetsUtil::RunTauJets(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchTauJet,const vector<fastjet::PseudoJet> & sorted_jets,const vector<PhysicsTower> & towers, const vector<TLorentzVector> & TrackCentral)
{
  TRootTauJet *elementTauJet;
  TLorentzVector JET;
  for (unsigned int i = 0; i < sorted_jets.size(); i++) {
    JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
    // Tau jet identification : 1! track and electromagnetic collimation
    if(fabs(JET.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) {
      double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
      if(
	 ( Energie_tau_central/JET.E() > DET->TAU_energy_frac ) &&
	 ( DET->NumTracks(TrackCentral,DET->TAU_track_pt,JET.Eta(),JET.Phi()) == 1 ) &&
	 ( JET.Pt() > DET->PTCUT_taujet)
	 ) {
	elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
	elementTauJet->Set(JET);
      } // if tau jet
    } // if JET.eta < tracker - tau_cone : Tau jet identification
    
  } // for itJet : loop on all jets
  
  
}