AliEmcalJet.cxx

Go to the documentation of this file.

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
#include "AliEmcalJet.h"

#include "AliLog.h"
#include "Riostream.h"

#include "AliParticleContainer.h"
#include "AliClusterContainer.h"

ClassImp(AliEmcalJet);

AliEmcalJet::AliEmcalJet() :
  AliVParticle(),
  fPt(0),
  fEta(0),
  fPhi(0),
  fM(0),
  fNEF(0),
  fArea(0),
  fAreaEta(0),
  fAreaPhi(0),
  fAreaE(0),
  fAreaEmc(-1),
  fAxisInEmcal(0),
  fFlavourTagging(0),
  fFlavourTracks(0),
  fMaxCPt(0),
  fMaxNPt(0),
  fMCPt(0),
  fNn(0),
  fNch(0),
  fPtEmc(0),
  fNEmc(0),
  fClusterIDs(),
  fTrackIDs(),
  fMatched(2),
  fMatchingType(0),
  fTaggedJet(0x0),
  fTagStatus(-1),
  fPtSub(0),
  fPtSubVect(0),
  fTriggers(0),
  fLabel(-1),
  fHasGhost(kFALSE),
  fGhosts(),
  fJetShapeProperties(0),
  fJetAcceptanceType(0),
  fParticleConstituents(),
  fClusterConstituents()
{
  fClosestJets[0] = 0;
  fClosestJets[1] = 0;
  fClosestJetsDist[0] = 999;
  fClosestJetsDist[1] = 999;
}

AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
  AliVParticle(),
  fPt(TMath::Sqrt(px * px + py* py)),
  fEta(TMath::ASinH(pz / fPt)),
  fPhi(0),
  fM(0),
  fNEF(0),
  fArea(0),
  fAreaEta(0),
  fAreaPhi(0),
  fAreaE(0),
  fAreaEmc(-1),
  fAxisInEmcal(0),
  fFlavourTagging(0),
  fFlavourTracks(0),
  fMaxCPt(0),
  fMaxNPt(0),
  fMCPt(0),
  fNn(0),
  fNch(0),
  fPtEmc(0),
  fNEmc(0),
  fClusterIDs(),
  fTrackIDs(),
  fMatched(2),
  fMatchingType(0),
  fTaggedJet(0x0),
  fTagStatus(-1),
  fPtSub(0),
  fPtSubVect(0),
  fTriggers(0),
  fLabel(-1),
  fHasGhost(kFALSE),
  fGhosts(),
  fJetShapeProperties(0),
  fJetAcceptanceType(0),
  fParticleConstituents(),
  fClusterConstituents()
{
  if (fPt != 0) {
    fPhi = TVector2::Phi_0_2pi(TMath::ATan2(py, px));
  }

  fClosestJets[0] = 0;
  fClosestJets[1] = 0;
  fClosestJetsDist[0] = 999;
  fClosestJetsDist[1] = 999;
}

AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
  AliVParticle(),
  fPt(pt),
  fEta(eta),
  fPhi(phi),
  fM(m),
  fNEF(0),
  fArea(0),
  fAreaEta(0),
  fAreaPhi(0),
  fAreaE(0),
  fAreaEmc(-1),
  fAxisInEmcal(0),
  fFlavourTagging(0),
  fFlavourTracks(0),
  fMaxCPt(0),
  fMaxNPt(0),
  fMCPt(0),
  fNn(0),
  fNch(0),
  fPtEmc(0),
  fNEmc(0),
  fClusterIDs(),
  fTrackIDs(),
  fMatched(2),
  fMatchingType(0),
  fTaggedJet(0x0),
  fTagStatus(-1),
  fPtSub(0),
  fPtSubVect(0),
  fTriggers(0),
  fLabel(-1),
  fHasGhost(kFALSE),
  fGhosts(),
  fJetShapeProperties(0),
  fJetAcceptanceType(0),
  fParticleConstituents(),
  fClusterConstituents()
{
  fPhi = TVector2::Phi_0_2pi(fPhi);

  fClosestJets[0] = 0;
  fClosestJets[1] = 0;
  fClosestJetsDist[0] = 999;
  fClosestJetsDist[1] = 999;
}

AliEmcalJet::AliEmcalJet(const AliEmcalJet& jet) :
  AliVParticle(jet),
  fPt(jet.fPt),
  fEta(jet.fEta),
  fPhi(jet.fPhi),
  fM(jet.fM),
  fNEF(jet.fNEF),
  fArea(jet.fArea),
  fAreaEta(jet.fAreaEta),
  fAreaPhi(jet.fAreaPhi),
  fAreaE(jet.fAreaE),
  fAreaEmc(jet.fAreaEmc),
  fAxisInEmcal(jet.fAxisInEmcal),
  fFlavourTagging(jet.fFlavourTagging),
  fFlavourTracks((jet.fFlavourTracks) ? new TObjArray(*(jet.fFlavourTracks)) : 0),
  fMaxCPt(jet.fMaxCPt),
  fMaxNPt(jet.fMaxNPt),
  fMCPt(jet.fMCPt),
  fNn(jet.fNn),
  fNch(jet.fNch),
  fPtEmc(jet.fPtEmc),
  fNEmc(jet.fNEmc),
  fClusterIDs(jet.fClusterIDs),
  fTrackIDs(jet.fTrackIDs),
  fMatched(jet.fMatched),
  fMatchingType(jet.fMatchingType),
  fTaggedJet(jet.fTaggedJet),
  fTagStatus(jet.fTagStatus),
  fPtSub(jet.fPtSub),
  fPtSubVect(jet.fPtSubVect),
  fTriggers(jet.fTriggers),
  fLabel(jet.fLabel),
  fHasGhost(jet.fHasGhost),
  fGhosts(jet.fGhosts),
  fJetShapeProperties(0),
  fJetAcceptanceType(jet.fJetAcceptanceType),
  fParticleConstituents(jet.fParticleConstituents),
  fClusterConstituents(jet.fClusterConstituents)

{
  // Copy constructor.
  fClosestJets[0]     = jet.fClosestJets[0];
  fClosestJets[1]     = jet.fClosestJets[1];
  fClosestJetsDist[0] = jet.fClosestJetsDist[0];
  fClosestJetsDist[1] = jet.fClosestJetsDist[1];

  if (jet.fJetShapeProperties) {
    fJetShapeProperties = new AliEmcalJetShapeProperties(*(jet.fJetShapeProperties));
  }
}

AliEmcalJet::~AliEmcalJet()
{
  if (fJetShapeProperties) delete fJetShapeProperties;
}

AliEmcalJet& AliEmcalJet::operator=(const AliEmcalJet& jet)
{
  if (this != &jet) {
    AliVParticle::operator=(jet);
    fPt                 = jet.fPt;
    fEta                = jet.fEta;
    fPhi                = jet.fPhi;
    fM                  = jet.fM;
    fNEF                = jet.fNEF;
    fArea               = jet.fArea;
    fAreaEta            = jet.fAreaEta;
    fAreaPhi            = jet.fAreaPhi;
    fAreaE              = jet.fAreaE;
    fAreaEmc            = jet.fAreaEmc;
    fAxisInEmcal        = jet.fAxisInEmcal;
    fFlavourTagging     = jet.fFlavourTagging;
    fFlavourTracks      = (jet.fFlavourTracks) ? new TObjArray(*(jet.fFlavourTracks)) : 0;
    fMaxCPt             = jet.fMaxCPt;
    fMaxNPt             = jet.fMaxNPt;
    fMCPt               = jet.fMCPt;
    fNn                 = jet.fNn;
    fNch                = jet.fNch;
    fPtEmc              = jet.fPtEmc;
    fNEmc               = jet.fNEmc;
    fClusterIDs         = jet.fClusterIDs;
    fTrackIDs           = jet.fTrackIDs;
    fClosestJets[0]     = jet.fClosestJets[0];
    fClosestJets[1]     = jet.fClosestJets[1];
    fClosestJetsDist[0] = jet.fClosestJetsDist[0];
    fClosestJetsDist[1] = jet.fClosestJetsDist[1];
    fMatched            = jet.fMatched;
    fTaggedJet          = jet.fTaggedJet;
    fTagStatus          = jet.fTagStatus;
    fPtSub              = jet.fPtSub;
    fPtSubVect          = jet.fPtSubVect;
    fTriggers           = jet.fTriggers;
    fLabel              = jet.fLabel;
    fHasGhost = jet.fHasGhost;
    fGhosts   = jet.fGhosts;
    if (jet.fJetShapeProperties) {
      fJetShapeProperties = new AliEmcalJetShapeProperties(*(jet.fJetShapeProperties));
    }
    fJetAcceptanceType  = jet.fJetAcceptanceType;
    fParticleConstituents = jet.fParticleConstituents;
    fClusterConstituents = jet.fClusterConstituents;
  }

  return *this;
}

Int_t AliEmcalJet::Compare(const TObject* obj) const
{
  //Return -1 if this is smaller than obj, 0 if objects are equal and 1 if this is larger than obj.

  if (obj == this) return 0;

  const AliEmcalJet* jet = dynamic_cast<const AliEmcalJet*>(obj);
  if (!jet) return 0;

  if (Pt() > jet->Pt()) return -1;
  else if (Pt() < jet->Pt()) return 1;
  else return 0;
}

void AliEmcalJet::GetMomentum(TLorentzVector& vec) const
{
  vec.SetPtEtaPhiE(fPt, fEta, fPhi, E());
}

Double_t AliEmcalJet::PtSub(Double_t rho, Bool_t save)
{
  Double_t ptcorr = fPt - rho * fArea;
  if (save) fPtSub = ptcorr;
  return ptcorr;
}

Double_t AliEmcalJet::PtSubVect(Double_t rho, Bool_t save)
{
  Double_t dx = Px() - rho * fArea * TMath::Cos(fAreaPhi);
  Double_t dy = Py() - rho * fArea * TMath::Sin(fAreaPhi);
  //Double_t dz = Pz() - rho * fArea * TMath::SinH(fAreaEta);
  Double_t ptcorr = TMath::Sqrt(dx * dx + dy * dy);
  if (save) fPtSubVect = ptcorr;
  return ptcorr;
}

TLorentzVector AliEmcalJet::SubtractRhoVect(Double_t rho, Bool_t save)
{
  TLorentzVector vecCorr;
  GetMomentum(vecCorr);
  TLorentzVector vecBg;
  vecBg.SetPtEtaPhiE(fArea, fAreaEta, fAreaPhi, fAreaE);
  vecBg *= rho;
  vecCorr -= vecBg;
  if (save) {
    Double_t dPhi = TMath::Abs(TVector2::Phi_mpi_pi(Phi() - vecCorr.Phi()));
    Int_t signum = dPhi <= TMath::PiOver2() ? 1 : -1;
    fPtSubVect = signum * vecCorr.Pt();
  }
  return vecCorr;
}

void AliEmcalJet::SortConstituents()
{
  std::sort(fClusterIDs.GetArray(), fClusterIDs.GetArray() + fClusterIDs.GetSize());
  std::sort(fTrackIDs.GetArray(), fTrackIDs.GetArray() + fTrackIDs.GetSize());
}

Double_t AliEmcalJet::DeltaR(const AliVParticle* part) const
{
  Double_t dPhi = Phi() - part->Phi();
  Double_t dEta = Eta() - part->Eta();
  dPhi = TVector2::Phi_mpi_pi(dPhi);
  return TMath::Sqrt(dPhi * dPhi + dEta * dEta);
}

std::vector<int> AliEmcalJet::GetPtSortedTrackConstituentIndexes(TClonesArray* tracks) const
{
  typedef std::pair<Double_t, Int_t> ptidx_pair;

  // Create vector for Pt sorting
  std::vector<ptidx_pair> pair_list;

  for (Int_t i_entry = 0; i_entry < GetNumberOfTracks(); i_entry++) {
    AliVParticle* track = Track(i_entry);
    if (!track) {
      AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)", i_entry, tracks->GetName(), TrackAt(i_entry), tracks->GetEntriesFast()));
      continue;
    }
    pair_list.push_back(std::make_pair(track->Pt(), i_entry));
  }

  std::stable_sort(pair_list.begin() , pair_list.end() , sort_descend());

  // return a vector of indexes of constituents (sorted descending by pt)
  std::vector <int> index_sorted_list;

  // populating the return object with indexes of sorted tracks
  for (auto it : pair_list) index_sorted_list.push_back(it.second);

  return index_sorted_list;
}

Double_t AliEmcalJet::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz) const
{
  Double_t pJetSq = P();
  pJetSq *= pJetSq;

  if (pJetSq > 1e-6) {
    return (trkPx * Px() + trkPy * Py() + trkPz * Pz()) / pJetSq;
  }
  else {
    AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f", GetName(), pJetSq));
    return -1;
  }
}

Double_t AliEmcalJet::GetZ(const AliVParticle* trk) const
{
  return GetZ(trk->Px(), trk->Py(), trk->Pz());
}

AliVParticle* AliEmcalJet::GetLeadingTrack(TClonesArray* tracks) const
{
  AliVParticle* maxTrack = 0;
  for (Int_t i = 0; i < GetNumberOfTracks(); i++) {
    AliVParticle* track = Track(i);
    if (!track) {
      AliError(Form("Unable to find jet track %d (global index %d) in the jet", i, TrackAt(i)));
      continue;
    }
    if (!maxTrack || track->Pt() > maxTrack->Pt())
      maxTrack = track;
  }

  return maxTrack;
}

AliVCluster* AliEmcalJet::GetLeadingCluster(TClonesArray* clusters) const
{
  AliVCluster* maxCluster = 0;
  for (Int_t i = 0; i < GetNumberOfClusters(); i++) {
    AliVCluster* cluster = Cluster(i);
    if (!cluster) {
      AliError(Form("Unable to find jet cluster %d (global index %d) in the jet", i, ClusterAt(i)));
      continue;
    }
    if (!maxCluster || cluster->E() > maxCluster->E())
      maxCluster = cluster;
  }

  return maxCluster;
}

void AliEmcalJet::ResetMatching()
{
  fClosestJets[0] = 0;
  fClosestJets[1] = 0;
  fClosestJetsDist[0] = 999;
  fClosestJetsDist[1] = 999;
  fMatched = 2;
}

Int_t AliEmcalJet::ContainsTrack(Int_t it) const
{
  for (Int_t i = 0; i < fTrackIDs.GetSize(); i++) {
    if (it == fTrackIDs[i]) return i;
  }
  return -1;
}

Int_t AliEmcalJet::ContainsCluster(Int_t ic) const
{
  for (Int_t i = 0; i < fClusterIDs.GetSize(); i++) {
    if (ic == fClusterIDs[i]) return i;
  }
  return -1;
}

TString AliEmcalJet::toString() const {
  return TString::Format("Jet pT = %.2f, eta = %.2f, phi = %.2f, max charged pT = %.2f, max neutral pT = %.2f, N tracks = %d, N clusters = %d, Area = %.2f, NEF = %.2f",
         Pt(), Eta(), Phi(), MaxChargedPt(), MaxNeutralPt(), GetNumberOfTracks(), GetNumberOfClusters(), Area(), NEF());
}

void AliEmcalJet::Print(Option_t* /*opt*/) const
{
  Printf("%s\n", toString().Data());
}

std::ostream &AliEmcalJet::Print(std::ostream &in) const {
  in << toString().Data();
  return in;
}

void AliEmcalJet::PrintConstituents(TClonesArray* tracks, TClonesArray* clusters) const
{
  for (Int_t i = 0; i < GetNumberOfTracks(); i++) {
    AliVParticle* part = Track(i);
    if (part) {
      Printf("Track %d (index = %d) pT = %.2f, eta = %.2f, phi = %.2f, PDG code = %d", i, TrackAt(i), part->Pt(), part->Eta(), part->Phi(), part->PdgCode());
    }
  }

  for (Int_t i = 0; i < GetNumberOfClusters(); i++) {
    AliVCluster* clus = Cluster(i);
    if (clus) {
      Printf("Cluster %d (index = %d) E = %.2f", i, ClusterAt(i), clus->E());
    }
  }
}

Double_t AliEmcalJet::GetFlavourTrackZ(Int_t i)  const
{
  if (P() < 1e-6) return 0.;
  AliVParticle* hftrack = GetFlavourTrack(i);
  return hftrack != 0 ? hftrack->P() / P() : 0.;
}

std::ostream &operator<<(std::ostream &in, const AliEmcalJet &myjet) {
  std::ostream &result = myjet.Print(in);
  return result;
}

void AliEmcalJet::AddGhost(const Double_t dPx, const Double_t dPy, const Double_t dPz, const Double_t dE)
{
  TLorentzVector ghost(dPx, dPy, dPz, dE);
  fGhosts.push_back(ghost);
  if (!fHasGhost) fHasGhost = kTRUE;
  return;
}

void AliEmcalJet::Clear(Option_t */*option*/)
{
  fClusterIDs.Set(0);
  fTrackIDs.Set(0);
  fClosestJets[0] = 0;
  fClosestJets[1] = 0;
  fClosestJetsDist[0] = 0;
  fClosestJetsDist[1] = 0;
  fMatched = 0;
  fPtSub = 0;
  fGhosts.clear();
  fHasGhost = kFALSE;
  fClusterConstituents.clear();
  fParticleConstituents.clear();
}

AliVParticle* AliEmcalJet::Track(Int_t idx) const
{
  return AliParticleContainer::GetEmcalContainerIndexMap().GetObjectFromGlobalIndex(TrackAt(idx));
}

AliVParticle* AliEmcalJet::TrackAt(Int_t idx, TClonesArray *ta) const
{
  if (!ta) return 0;
  auto res =  AliParticleContainer::GetEmcalContainerIndexMap().LocalIndexFromGlobalIndex(TrackAt(idx));
  if (res.second != ta) {
    AliWarning(Form("TClonesArray %s that was passed does not correspond to the passed index! The index belongs to a different TClonesArray named %s! Returning the object corresponding to the index (not the passed TClonesArray)! Consider fixing by updating to jet->Track(index).", ta->GetName(), res.second->GetName()));
  }
  return dynamic_cast<AliVParticle*>(res.second->At(res.first));
}

Int_t AliEmcalJet::ContainsTrack(AliVParticle* track, TClonesArray* tracks) const
{
  if (!tracks || !track) return 0;
  return ContainsTrack(tracks->IndexOf(track));
}

AliVCluster* AliEmcalJet::Cluster(Int_t idx) const
{
  return AliClusterContainer::GetEmcalContainerIndexMap().GetObjectFromGlobalIndex(ClusterAt(idx));
}

AliVCluster* AliEmcalJet::ClusterAt(Int_t idx, TClonesArray *ca) const
{
  if (!ca) return 0;

  auto res =  AliClusterContainer::GetEmcalContainerIndexMap().LocalIndexFromGlobalIndex(ClusterAt(idx));
  if (res.second != ca) {
    AliWarning(Form("TClonesArray %s that was passed does not correspond to the passed index! The index belongs to a different TClonesArray named %s! Returning the object corresponding to the index (not the passed TClonesArray)! Consider fixing by updating to jet->Cluster(index).", ca->GetName(), res.second->GetName()));
  }
  return dynamic_cast<AliVCluster*>(res.second->At(res.first));
}

Int_t AliEmcalJet::ContainsCluster(AliVCluster* cluster, TClonesArray* clusters) const
{
  if (!clusters || !cluster) return 0;
  return ContainsCluster(clusters->IndexOf(cluster));
}

const PWG::JETFW::AliEmcalClusterJetConstituent *AliEmcalJet::ClusterConstituentAt(unsigned int icl) const {
  return (icl < fClusterConstituents.size()) ? &fClusterConstituents[icl] : nullptr;
}

const PWG::JETFW::AliEmcalParticleJetConstituent *AliEmcalJet::ParticleConstituentAt(unsigned int ipart) const {
  return (ipart < fParticleConstituents.size()) ? &fParticleConstituents[ipart] : nullptr;
}

const PWG::JETFW::AliEmcalClusterJetConstituent *AliEmcalJet::GetLeadingClusterConstituent() const {
  if(!fClusterConstituents.size()) return nullptr;
  const PWG::JETFW::AliEmcalClusterJetConstituent *leading(nullptr);
  for(const auto &c : fClusterConstituents) {
    if(!leading) leading = &c;
    else {
      if(c.E() > leading->E()) leading = &c;
    }
  }
  return leading;
}

const PWG::JETFW::AliEmcalParticleJetConstituent *AliEmcalJet::GetLeadingParticleConstituent() const {
  if(!fParticleConstituents.size()) return nullptr;
  const PWG::JETFW::AliEmcalParticleJetConstituent *leading(nullptr);
  for(const auto &p : fParticleConstituents) {
    if(!leading) leading = &p;
    else {
      if(p.Pt() > leading->Pt()) leading = &p;
    }
  }
  return leading;
}

bool AliEmcalJet::HasClusterConstituent(const AliVCluster *const clust) const {
  PWG::JETFW::AliEmcalClusterJetConstituent test(clust);
  bool found(false);
  for(auto c : fClusterConstituents) {
    if(test == c) {
      found = true;
      break;
    }
  }
  return found;
}

bool AliEmcalJet::HasParticleConstituent(const AliVParticle *const part) const {
  PWG::JETFW::AliEmcalParticleJetConstituent test(part);
  bool found(false);
  for(auto p : fParticleConstituents) {
    if(test == p) {
      found = true;
      break;
    }
  }
  return found;
}


Double_t AliEmcalJet::GetXi(const AliVParticle* trk) const
{
  return TMath::Log(1 / GetZ(trk));
}

Double_t AliEmcalJet::GetXi( const Double_t trkPx, const Double_t trkPy, const Double_t trkPz ) const
{
  return TMath::Log(1 / GetZ(trkPx, trkPy, trkPz));
}

void AliEmcalJet::AddFlavourTrack(AliVParticle* hftrack)
{
  if (!fFlavourTracks) fFlavourTracks = new TObjArray();
  fFlavourTracks->Add(hftrack);
}

void AliEmcalJet::AddParticleConstituent(const AliVParticle *const part, Bool_t isEmbedding, UInt_t globalIndex) {
  PWG::JETFW::AliEmcalParticleJetConstituent constituent(part);
  constituent.SetIsFromEmbeddedEvent(isEmbedding);
  constituent.SetGlobalIndex(globalIndex);
  AddParticleConstituent(constituent);
}

void AliEmcalJet::AddParticleConstituent(const PWG::JETFW::AliEmcalParticleJetConstituent &part){
  fParticleConstituents.emplace_back(part);
}


void AliEmcalJet::AddClusterConstituent(const AliVCluster *const clust, AliVCluster::VCluUserDefEnergy_t endef, Double_t *pvec, Bool_t isEmbedding, UInt_t globalIndex) {
  PWG::JETFW::AliEmcalClusterJetConstituent constituent(clust, endef, pvec);
  constituent.SetIsFromEmbeddedEvent(isEmbedding);
  constituent.SetGlobalIndex(globalIndex);
  AddClusterConstituent(constituent);
}

void AliEmcalJet::AddClusterConstituent(const PWG::JETFW::AliEmcalClusterJetConstituent &clust) {
  fClusterConstituents.emplace_back(clust);
}

AliVParticle* AliEmcalJet::GetFlavourTrack(Int_t i) const
{
  if (!fFlavourTracks || i < 0 || i >= fFlavourTracks->GetEntriesFast()) return 0;

  return  static_cast<AliVParticle*>(fFlavourTracks->At(i));
}

AliVParticle* AliEmcalJet::RemoveFlavourTrack(Int_t i)
{
  if (!fFlavourTracks || i < 0 || i >= fFlavourTracks->GetEntriesFast()) return 0;

  return static_cast<AliVParticle*>(fFlavourTracks->RemoveAt(i));
}