Index: /trunk/interface/JetsUtil.h
===================================================================
--- /trunk/interface/JetsUtil.h	(revision 215)
+++ /trunk/interface/JetsUtil.h	(revision 215)
@@ -0,0 +1,75 @@
+#ifndef _JETSUTIL_H_
+#define _JETSUTIL_H_
+
+/*
+  ---- 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
+*/
+
+/// \file SmearUtil.h
+/// \brief RESOLution class, and some generic definitions
+
+
+#include <vector>
+#include "TLorentzVector.h"
+
+#include "SmearUtil.h"
+#include "BlockClasses.h"
+#include "TSimpleArray.h"
+#include "ExRootTreeReader.h"
+#include "ExRootTreeWriter.h"
+#include "ExRootTreeBranch.h"
+#include "PhysicsTower.hh"
+#include "TRandom.h"
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <iomanip>
+#include <cstring>
+#include "PseudoJet.hh"
+#include "ClusterSequence.hh"
+
+// get info on how fastjet was configured
+#include "Utilities/Fastjet/include/fastjet/config.h"
+
+// make sure we have what is needed
+#ifdef ENABLE_PLUGIN_SISCONE
+#  include "SISConePlugin.hh"
+#endif
+#ifdef ENABLE_PLUGIN_CDFCONES
+#  include "CDFMidPointPlugin.hh"
+#  include "CDFJetCluPlugin.hh"
+#endif
+
+using namespace std;
+
+class JetsUtil
+{
+ public:
+  /// Constructor
+  JetsUtil();
+  JetsUtil(const string& DetDatacard);
+  JetsUtil(const RESOLution* DetDatacard);
+  JetsUtil(const JetsUtil& jet);
+  JetsUtil& operator=(const JetsUtil& jet);
+  ~JetsUtil() {delete DET; delete plugins;};
+  void init(); // for constructors
+
+  vector<fastjet::PseudoJet> RunJets(const vector<fastjet::PseudoJet>&  input_particles);
+  void RunJetBtagging(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchJet,const vector<fastjet::PseudoJet> & sorted_jets,const TSimpleArray<TRootGenParticle> & NFCentralQ);
+  void RunTauJets(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchTauJet,const vector<fastjet::PseudoJet> & sorted_jets,const vector<PhysicsTower> & towers, const vector<TLorentzVector> & TrackCentral);
+
+ private:
+  RESOLution *DET;
+  fastjet::JetDefinition jet_def;
+  fastjet::JetDefinition::Plugin * plugins;
+  vector<fastjet::PseudoJet> inclusive_jets;
+  vector<fastjet::PseudoJet> sorted_jets;
+};
+
+#endif
Index: /trunk/src/JetsUtil.cc
===================================================================
--- /trunk/src/JetsUtil.cc	(revision 215)
+++ /trunk/src/JetsUtil.cc	(revision 215)
@@ -0,0 +1,163 @@
+/*
+ *   ---- 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
+  
+  
+}
+
