source: svn/trunk/src/JetsUtil.cc@ 814

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  Delphes, a framework for the fast simulation  |         **
**         |       of a  generic collider experiment        |         **
**          \-------------  arXiv:0903.2225v1  ------------/          **
**                                                                    **
**                                                                    **
**   This package uses:                                               **
**   ------------------                                               **
**   ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86        **
**   FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]      **
**   Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]       **
**   FROG: [hep-ex/0901.2718v1]                                       **
**   HepMC: Comput. Phys. Commun.134 (2001) 41                        **
**                                                                    **
** ------------------------------------------------------------------ **
**                                                                    **
**   Main authors:                                                    **
**   -------------                                                    **
**                                                                    **
**                Severine Ovyn                Xavier Rouby           **
**          severine.ovyn@uclouvain.be      xavier.rouby@cern         **
**                                                                    **
**         Center for Particle Physics and Phenomenology (CP3)        **
**                Universite catholique de Louvain (UCL)              **
**                     Louvain-la-Neuve, Belgium                      **
**                                                                    **
**                      Copyright (C) 2008-2009,                      **
**                        All rights reserved.                        **
**                                                                    **
***********************************************************************/

#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,  const vector<TRootTracks> & TrackCentral,  vector<int> &NTrackJet, vector<float> &EHADEEM,D_CaloTowerList list_of_active_towers, vector<int> &NCALO)
{
  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);
      //Bin tracks to make the link 
      float  iEtaTrack[TrackCentral.size()];
      float  iPhiTrack[TrackCentral.size()];
      for(unsigned int t = 0; t < TrackCentral.size(); t++)
        {
          if(!DET->JET_Eflow) 
            DET->BinEtaPhi(TrackCentral[t].PhiOuter,TrackCentral[t].EtaOuter,iPhiTrack[t],iEtaTrack[t]); 
          else {
            iPhiTrack[t] = TrackCentral[t].Phi;    iEtaTrack[t] = TrackCentral[t].Eta;  // should NOT be PhiOuter and EtaOuter
          }
        }
      int numTrackJet;
      for (unsigned int i = 0; i < sorted_jets.size(); i++) 
          {
            //check number of tracks in jets 
            numTrackJet=0;
            vector<fastjet::PseudoJet> constituents = clust_seq.constituents(sorted_jets[i]);
            for(unsigned int it = 0; it < TrackCentral.size(); it++)
               {
                 for (unsigned int j = 0; j < constituents.size(); j++)
                     {
                       if(DeltaR(iPhiTrack[it], iEtaTrack[it], constituents[j].phi(), constituents[j].eta())<0.001)numTrackJet++;
                     }
               }
            NTrackJet.push_back(numTrackJet);
            //now, get EHad over EEm
            float EmVal=0,HadVal=0;
            NCALO.push_back(constituents.size());
            for (unsigned int j = 0; j < constituents.size(); j++)
                {
                  //D_CaloTower calConsti(list_of_active_towers.getElement(constituents[j].eta(),constituents[j].phi()));
                  // bug fix! bad association of phi [-pi;pi] <-> [0 ; 2pi]
                  D_CaloTower calConsti(list_of_active_towers.getElement(constituents[j].eta(),constituents[j].phi_std()));
                  EmVal  += calConsti.getEem();
                  HadVal += calConsti.getEhad();
                }
            EHADEEM.push_back(HadVal/EmVal);
       }
    }
  
  return sorted_jets;                                                    
}

vector<fastjet::PseudoJet> JetsUtil::RunJetsResol(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<TRootC::GenParticle>& NFCentralQ, const vector<int> &NTrackJet, const vector<float> &EHADEEM,  const vector<int> &NCALO)
{
  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);
        elementJet->NTracks = NTrackJet[i];
        elementJet->EHoverEE = EHADEEM[i];
        elementJet->NCalo = NCALO[i];

        // 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<TRootTracks> & TrackCentral, const vector<int> &NTrackJet, const vector<float> &EHADEEM,const vector<int> &NCALO)
{
  TRootTauJet *elementTauJet;
  TLorentzVector JET;
  float charge=0;

  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! or 3! 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());
      int NumTrackTau = DET->NumTracks(charge,TrackCentral,DET->TAU_track_pt,JET.Eta(),JET.Phi());
      if(
         ( Energie_tau_central/JET.E() > DET->TAU_energy_frac ) &&
         (( NumTrackTau == 1 ) || ( NumTrackTau == 3 )) &&
         ( JET.Pt() > DET->PTCUT_taujet)
         ) {
        elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
        elementTauJet->Set(JET);
        elementTauJet->NTracks=NTrackJet[i];
        elementTauJet->Charge = charge;
        elementTauJet->EHoverEE = EHADEEM[i];
        elementTauJet->NCalo = NCALO[i];
      } // if tau jet
    } // if JET.eta < tracker - tau_cone : Tau jet identification
  } // for itJet : loop on all jets
}