/*
 *   ---- 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 "interface/JetUtils.h"

using namespace std;

JetsUtil::JetsUtil() {
  
  switch(JETALGO) {
  default:
  case 1: {
    
    // set up a CDF midpoint jet definition
#ifdef ENABLE_PLUGIN_CDFCONES
    plugins = new fastjet::CDFJetCluPlugin(SEEDTHRESHOLD,CONERADIUS,C_ADJACENCYCUT,C_MAXITERATIONS,C_IRATCH,OVERLAPTHRESHOLD);
    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 (SEEDTHRESHOLD,CONERADIUS,M_CONEAREAFRACTION,M_MAXPAIRSIZE,M_MAXITERATIONS,OVERLAPTHRESHOLD);
    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 (CONERADIUS,OVERLAPTHRESHOLD,NPASS, PROTOJET_PTMIN);
    jet_def = fastjet::JetDefinition(plugins);
#else
    plugins = NULL;
#endif
  }
    break;
    
  case 4: {
    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, CONERADIUS);
  }
    break;
    
  case 5: {
    jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,CONERADIUS);
  }
    break;
    
  case 6: {
    jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,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() > JET_pt)
      {
	elementJet = (TRootJet*) branchJet->NewEntry();
	elementJet->Set(JET);
	// b-jets
	bool btag=false;
	if((fabs(JET.Eta()) < MAX_TRACKER && 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()) < (MAX_TRACKER - TAU_CONE_TRACKS)) {
      double Energie_tau_central = EnergySmallCone(towers,JET.Eta(),JET.Phi());
      if(
	 ( Energie_tau_central/JET.E() > TAU_EM_COLLIMATION ) &&
	 ( NumTracks(TrackCentral,PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 ) &&
	 ( JET.Pt() > TAUJET_pt)
	 ) {
	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
  
  
}