[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliEmcalJetTask.cxx

//
// Emcal jet finder task.
//
// Authors: C.Loizides, S.Aiola

#include <vector>
#include "AliEmcalJetTask.h"

#include <TChain.h>
#include <TClonesArray.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TMath.h>
#include <TRandom3.h>

#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliEMCALGeometry.h"
#include "AliESDEvent.h"
#include "AliEmcalJet.h"
#include "AliEmcalParticle.h"
#include "AliFJWrapper.h"
#include "AliMCEvent.h"
#include "AliVCluster.h"
#include "AliVEvent.h"
#include "AliVParticle.h"
using std::cout;
using std::endl;
using std::cerr;

ClassImp(AliEmcalJetTask)

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask() : 
  AliAnalysisTaskSE("AliEmcalJetTask"),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fJetType(kNone),
  fConstSel(0),
  fMCConstSel(0),
  fMarkConst(kFALSE),
  fRadius(0.4),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fPhiMin(0),
  fPhiMax(TMath::TwoPi()),
  fEtaMin(-0.9),
  fEtaMax(+0.9),
  fMinJetArea(0.005),
  fMinJetPt(1.0),
  fJetPhiMin(-10),
  fJetPhiMax(+10),
  fJetEtaMin(-1),
  fJetEtaMax(+1),
  fGhostArea(0.005),
  fMinMCLabel(0),
  fRecombScheme(fastjet::pt_scheme),
  fTrackEfficiency(1.),
  fIsInit(0),
  fIsPSelSet(0),
  fIsMcPart(0),
  fIsEmcPart(0),
  fLegacyMode(kFALSE),
  fCodeDebug(kFALSE),
  fJets(0),
  fEvent(0),
  fTracks(0),
  fClus(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask(const char *name) : 
  AliAnalysisTaskSE(name),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fJetType(kAKT|kFullJet|kRX1Jet),
  fConstSel(0),
  fMCConstSel(0),
  fMarkConst(kFALSE),
  fRadius(0.4),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fPhiMin(0),
  fPhiMax(TMath::TwoPi()),
  fEtaMin(-0.9),
  fEtaMax(+0.9),
  fMinJetArea(0.001),
  fMinJetPt(1.0),
  fJetPhiMin(-10),
  fJetPhiMax(+10),
  fJetEtaMin(-1),
  fJetEtaMax(+1),
  fGhostArea(0.005),
  fMinMCLabel(0),
  fRecombScheme(fastjet::pt_scheme),
  fTrackEfficiency(1.),
  fIsInit(0),
  fIsPSelSet(0),
  fIsMcPart(0),
  fIsEmcPart(0),
  fLegacyMode(kFALSE),
  fCodeDebug(kFALSE),
  fJets(0),
  fEvent(0),
  fTracks(0),
  fClus(0)
{
  // Standard constructor.

  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
}

//________________________________________________________________________
AliEmcalJetTask::~AliEmcalJetTask()
{
  // Destructor
}

//________________________________________________________________________
void AliEmcalJetTask::UserCreateOutputObjects()
{
  // Create user objects.

  fJets = new TClonesArray("AliEmcalJet");
  fJets->SetName(fJetsName);
}

//________________________________________________________________________
void AliEmcalJetTask::UserExec(Option_t *) 
{
  // Main loop, called for each event.
  if (!fIsInit) {
    if (!DoInit())
      return;
    fIsInit = kTRUE;
  }

  // clear the jet array (normally a null operation)
  fJets->Delete();

  FindJets();
}

//________________________________________________________________________
void AliEmcalJetTask::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}

//________________________________________________________________________
void AliEmcalJetTask::FindJets()
{
  // Find jets.
  if (!fTracks && !fClus){
    cout << "WARNING NO TRACKS OR CLUSTERS:"  <<endl;
    return;
  }

  TString name("kt");
  fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
  if ((fJetType & kAKT) != 0) {
    name  = "antikt";
    jalgo = fastjet::antikt_algorithm;
    AliDebug(1,"Using AKT algorithm");
  }
  else {
    AliDebug(1,"Using KT algorithm");
  }

  if ((fJetType & kR020Jet) != 0)
    fRadius = 0.2;
  else if ((fJetType & kR030Jet) != 0)
    fRadius = 0.3; 
  else if ((fJetType & kR040Jet) != 0)
    fRadius = 0.4;

  // setup fj wrapper
  AliFJWrapper fjw(name, name);
  fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
  fjw.SetGhostArea(fGhostArea);
  fjw.SetR(fRadius);
  fjw.SetAlgorithm(jalgo);
  fjw.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme)); 
  fjw.SetMaxRap(fEtaMax);
  fjw.Clear();

  // get primary vertex
  Double_t vertex[3] = {0, 0, 0};
  if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(vertex);

  AliDebug(2,Form("Jet type = %d", fJetType));

  if ((fIsMcPart || ((fJetType & kFullJet) != 0) || ((fJetType & kChargedJet) != 0)) && fTracks) {
    const Int_t Ntracks = fTracks->GetEntries();
    for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
      AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
      if (!t)
        continue;
      if (fIsMcPart) {
	if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
	  AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
	  continue;
	}
	if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
	  AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
	  continue;
	}
      }
      if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) {
	if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
	  AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
	  continue;
	}
	else {
	  AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
	}
      }
      else {
	if (t->TestBits(fConstSel) != (Int_t)fConstSel) {
	  AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
	  continue;
	}
	else {
	  AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));	    
	}
      }	    
      if (t->Pt() < fMinJetTrackPt) 
        continue;
      Double_t eta = t->Eta();
      Double_t phi = t->Phi();
      if ((eta<fEtaMin) || (eta>fEtaMax) ||
          (phi<fPhiMin) || (phi>fPhiMax))
        continue;

      // artificial inefficiency
      if (fTrackEfficiency < 1.) {
	Double_t rnd = gRandom->Rndm();
	if (fTrackEfficiency < rnd) {
	  AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", iTracks));
	  continue;
	}
      }

      // offset of 100 for consistency with cluster ids
      AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt()));
      fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
    }
  }

  if ((((fJetType & kFullJet) != 0) || ((fJetType & kNeutralJet) != 0)) && fClus) {
    const Int_t Nclus = fClus->GetEntries();
    for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
      AliVCluster *c = 0;
      Double_t cEta=0,cPhi=0,cPt=0;
      Double_t cPx=0,cPy=0,cPz=0;
      if (fIsEmcPart) {
	AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
	if (!ep)
	  continue;

	c = ep->GetCluster();
	if (!c)
	  continue;

	if (c->GetLabel() > fMinMCLabel) {
	  if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
	    AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
	    continue;
	  }
	  else {
	    AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
	  }
	}
	else {
	  if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
	    AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
	    continue;
	  }
	  else {
	    AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));	    
	  }
	}

	cEta = ep->Eta();
	cPhi = ep->Phi();
	cPt  = ep->Pt();
	cPx  = ep->Px();
	cPy  = ep->Py();
	cPz  = ep->Pz();
      } else {
	c = static_cast<AliVCluster*>(fClus->At(iClus));
	if (!c)
	  continue;

	if (c->GetLabel() > fMinMCLabel) {
	  if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
	    AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
	    continue;
	  }
	  else {
	    AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
	  }
	}
	else {
	  if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
	    AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
	    continue;
	  }
	  else {
	    AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));	    
	  }
	}

	TLorentzVector nP;
	c->GetMomentum(nP, vertex);
	cEta = nP.Eta();
	cPhi = nP.Phi();
	cPt  = nP.Pt();
	cPx  = nP.Px();
	cPy  = nP.Py();
	cPz  = nP.Pz();
      }
      if (!c->IsEMCAL())
	continue;
      if (cPt < fMinJetClusPt) 
	continue;
      if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
	  (cPhi<fPhiMin) || (cPhi>fPhiMax))
	continue;
      // offset of 100 to skip ghost particles uid = -1
      AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
      fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
    }
  }

  // setting legacy mode
  if (fLegacyMode) { 
    fjw.SetLegacyMode(kTRUE);
  }

  // run jet finder
  fjw.Run();

  // get geometry
  AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
  if (!geom) {
    AliFatal(Form("%s: Can not create geometry", GetName()));
    return;
  }

  // loop over fastjet jets
  std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
  // sort jets according to jet pt
  static Int_t indexes[9999] = {-1};
  GetSortedArray(indexes, jets_incl);

  AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
  for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
    Int_t ij = indexes[ijet];
    AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));

    if (jets_incl[ij].perp()<fMinJetPt) 
      continue;
    if (fjw.GetJetArea(ij)<fMinJetArea)
      continue;
    if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
	(jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
      continue;

    AliEmcalJet *jet = new ((*fJets)[jetCount]) 
      AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());

    // loop over constituents
    std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
    jet->SetNumberOfTracks(constituents.size());
    jet->SetNumberOfClusters(constituents.size());

    Int_t nt            = 0;
    Int_t nc            = 0;
    Double_t neutralE   = 0;
    Double_t maxCh      = 0;
    Double_t maxNe      = 0;
    Int_t gall          = 0;
    Int_t gemc          = 0;
    Int_t cemc          = 0;
    Int_t ncharged      = 0;
    Int_t nneutral      = 0;
    Double_t mcpt       = 0;
    Double_t emcpt      = 0;

    for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
      Int_t uid = constituents[ic].user_index();
      AliDebug(3,Form("Processing constituent %d", uid));
      if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
        ++gall;
        Double_t gphi = constituents[ic].phi();
        if (gphi<0) 
          gphi += TMath::TwoPi();
        gphi *= TMath::RadToDeg();
        Double_t geta = constituents[ic].eta();
        if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
            (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
          ++gemc; 
      }	else if ((uid > 0) && fTracks) { // track constituent
	Int_t tid = uid - 100;
        AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
        if (!t) {
	  AliError(Form("Could not find track %d",tid));
          continue;
	}
	if (jetCount < fMarkConst) {
	  AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
	  t->SetBit(fJetType);
	}
        Double_t cEta = t->Eta();
        Double_t cPhi = t->Phi();
        Double_t cPt  = t->Pt();
        Double_t cP   = t->P();
        if (t->Charge() == 0) {
          neutralE += cP;
          ++nneutral;
          if (cPt > maxNe)
            maxNe = cPt;
        } else {
          ++ncharged;
          if (cPt > maxCh)
            maxCh = cPt;
        }
        if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
          mcpt += cPt;

        if (cPhi<0) 
          cPhi += TMath::TwoPi();
        cPhi *= TMath::RadToDeg();
        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
          emcpt += cPt;
          ++cemc;
        }

        jet->AddTrackAt(tid, nt);
        ++nt;
      } else if (fClus) { // cluster constituent
	Int_t cid = -uid - 100;
	AliVCluster *c = 0;
        Double_t cEta=0,cPhi=0,cPt=0,cP=0;
        if (fIsEmcPart) {
          AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
          if (!ep)
            continue;
          c = ep->GetCluster();
          if (!c)
            continue;
	  if (jetCount < fMarkConst)
	    ep->SetBit(fJetType);
          cEta = ep->Eta();
          cPhi = ep->Phi();
          cPt  = ep->Pt();
          cP   = ep->P();
        } else {
          c = static_cast<AliVCluster*>(fClus->At(cid));
          if (!c)
            continue;
	  if (jetCount < fMarkConst)
	    c->SetBit(fJetType);
          TLorentzVector nP;
          c->GetMomentum(nP, vertex);
          cEta = nP.Eta();
          cPhi = nP.Phi();
          cPt  = nP.Pt();
          cP   = nP.P();
        }

        neutralE += cP;
        if (cPt > maxNe)
          maxNe = cPt;

        if (c->GetLabel() > fMinMCLabel) // MC particle
          mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;

        if (cPhi<0) 
          cPhi += TMath::TwoPi();
        cPhi *= TMath::RadToDeg();
        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
          emcpt += cPt;
          ++cemc;
        }

        jet->AddClusterAt(cid, nc);
        ++nc;
        ++nneutral;
      } else {
        AliError(Form("%s: No logical way to end up here.", GetName()));
        continue;
      }
    } /* end of constituent loop */

    jet->SetNumberOfTracks(nt);
    jet->SetNumberOfClusters(nc);
    jet->SortConstituents();
    jet->SetMaxNeutralPt(maxNe);
    jet->SetMaxChargedPt(maxCh);
    jet->SetNEF(neutralE / jet->E());
    fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
    jet->SetArea(area.perp());
    jet->SetAreaEta(area.eta());
    jet->SetAreaPhi(area.phi());
    jet->SetNumberOfCharged(ncharged);
    jet->SetNumberOfNeutrals(nneutral);
    jet->SetMCPt(mcpt);
    jet->SetNEmc(cemc);
    jet->SetPtEmc(emcpt);

    if (gall > 0)
      jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
    else 
      jet->SetAreaEmc(-1);
    if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && 
	(jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
	(jet->Eta() > geom->GetArm1EtaMin()) && 
	(jet->Eta() < geom->GetArm1EtaMax()))
      jet->SetAxisInEmcal(kTRUE);

    AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
    jetCount++;
  }
  //fJets->Sort();
}

//________________________________________________________________________
Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
{
  // Get the leading jets.

  static Float_t pt[9999] = {0};

  const Int_t n = (Int_t)array.size();

  if (n < 1)
    return kFALSE;
  
  for (Int_t i = 0; i < n; i++) 
    pt[i] = array[i].perp();

  TMath::Sort(n, pt, indexes);

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliEmcalJetTask::DoInit()
{
  // Init. Return true if successful.

  if (fTrackEfficiency < 1.) {
    if (gRandom) delete gRandom;
    gRandom = new TRandom3(0);
  }

  // get input collections
  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();

  // get the event
  fEvent = InputEvent();
  if (!fEvent) {
    AliError(Form("%s: Could not retrieve event! Returning", GetName()));
    return 0;
  }

  // add jets to event if not yet there
  if (!(fEvent->FindListObject(fJetsName)))
    fEvent->AddObject(fJets);
  else {
    AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
    return 0;
  }

  if (fTracksName == "Tracks")
    am->LoadBranch("Tracks");
  if (!fTracks && !fTracksName.IsNull()) {
    fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
    if (!fTracks) {
      AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
      return 0;
    }
  }
  if (fTracks) {
    TClass cls(fTracks->GetClass()->GetName());
    if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
      fIsMcPart = 1;
  }
  
  if (fCaloName == "CaloClusters")
    am->LoadBranch("CaloClusters");
  if (!fClus && !fCaloName.IsNull()) {
    fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
    if (!fClus) {
      AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
      return 0;
    }
  }
  if (fClus) {
    TClass cls(fClus->GetClass()->GetName());
    if (cls.InheritsFrom("AliEmcalParticle"))
      fIsEmcPart = 1;
  }

  return 1;
}
Commit	Line	Data
542ac102	1	//
	2	// Emcal jet finder task.
	3	//
b65fac7a	4	// Authors: C.Loizides, S.Aiola
48d874ff	5
c64cb1f6	6	#include <vector>
04905231	7	#include "AliEmcalJetTask.h"
04905231	8
6ea168d0	9	#include <TChain.h>
48d874ff	10	#include <TClonesArray.h>
6ea168d0	11	#include <TList.h>
6ea168d0	12	#include <TLorentzVector.h>
ca5c29fa	13	#include <TMath.h>
6421eeb0	14	#include <TRandom3.h>
b65fac7a	15
48d874ff	16	#include "AliAnalysisManager.h"
6ea168d0	17	#include "AliCentrality.h"
f5f3c8e9	18	#include "AliEMCALGeometry.h"
04905231	19	#include "AliESDEvent.h"
f5f3c8e9	20	#include "AliEmcalJet.h"
04905231	21	#include "AliEmcalParticle.h"
63206072	22	#include "AliFJWrapper.h"
04905231	23	#include "AliMCEvent.h"
35789a2d	24	#include "AliVCluster.h"
b65fac7a	25	#include "AliVEvent.h"
04905231	26	#include "AliVParticle.h"
c2aad3ae	27	using std::cout;
	28	using std::endl;
	29	using std::cerr;
48d874ff	30
914d486c	31	ClassImp(AliEmcalJetTask)
48d874ff	32
6ea168d0	33	//________________________________________________________________________
914d486c	34	AliEmcalJetTask::AliEmcalJetTask() :
914d486c	35	AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0	36	fTracksName("Tracks"),
	37	fCaloName("CaloClusters"),
	38	fJetsName("Jets"),
6a20534a	39	fJetType(kNone),
d1f2e0d9	40	fConstSel(0),
d1f2e0d9	41	fMCConstSel(0),
6a20534a	42	fMarkConst(kFALSE),
6ea168d0	43	fRadius(0.4),
914d486c	44	fMinJetTrackPt(0.15),
914d486c	45	fMinJetClusPt(0.15),
04905231	46	fPhiMin(0),
04905231	47	fPhiMax(TMath::TwoPi()),
83086e6c	48	fEtaMin(-0.9),
83086e6c	49	fEtaMax(+0.9),
cba7b46e	50	fMinJetArea(0.005),
f5f3c8e9	51	fMinJetPt(1.0),
c9ad772c	52	fJetPhiMin(-10),
83086e6c	53	fJetPhiMax(+10),
c9ad772c	54	fJetEtaMin(-1),
83086e6c	55	fJetEtaMax(+1),
cba7b46e	56	fGhostArea(0.005),
a7477843	57	fMinMCLabel(0),
751cb2dc	58	fRecombScheme(fastjet::pt_scheme),
6421eeb0	59	fTrackEfficiency(1.),
04905231	60	fIsInit(0),
d14d6e2c	61	fIsPSelSet(0),
04905231	62	fIsMcPart(0),
04905231	63	fIsEmcPart(0),
5d0a5007	64	fLegacyMode(kFALSE),
ffdca6aa	65	fCodeDebug(kFALSE),
b65fac7a	66	fJets(0),
04905231	67	fEvent(0),
	68	fTracks(0),
	69	fClus(0)
6ea168d0	70	{
6ea168d0	71	// Default constructor.
6ea168d0	72	}
6ea168d0	73
48d874ff	74	//________________________________________________________________________
914d486c	75	AliEmcalJetTask::AliEmcalJetTask(const char *name) :
63206072	76	AliAnalysisTaskSE(name),
7efbea04	77	fTracksName("Tracks"),
48d874ff	78	fCaloName("CaloClusters"),
7efbea04	79	fJetsName("Jets"),
6a20534a	80	fJetType(kAKT\|kFullJet\|kRX1Jet),
d1f2e0d9	81	fConstSel(0),
d1f2e0d9	82	fMCConstSel(0),
6a20534a	83	fMarkConst(kFALSE),
48d874ff	84	fRadius(0.4),
6ea168d0	85	fMinJetTrackPt(0.15),
914d486c	86	fMinJetClusPt(0.15),
04905231	87	fPhiMin(0),
04905231	88	fPhiMax(TMath::TwoPi()),
83086e6c	89	fEtaMin(-0.9),
83086e6c	90	fEtaMax(+0.9),
cba7b46e	91	fMinJetArea(0.001),
f5f3c8e9	92	fMinJetPt(1.0),
c9ad772c	93	fJetPhiMin(-10),
83086e6c	94	fJetPhiMax(+10),
c9ad772c	95	fJetEtaMin(-1),
83086e6c	96	fJetEtaMax(+1),
cba7b46e	97	fGhostArea(0.005),
a7477843	98	fMinMCLabel(0),
751cb2dc	99	fRecombScheme(fastjet::pt_scheme),
6421eeb0	100	fTrackEfficiency(1.),
04905231	101	fIsInit(0),
d14d6e2c	102	fIsPSelSet(0),
04905231	103	fIsMcPart(0),
04905231	104	fIsEmcPart(0),
31bab07e	105	fLegacyMode(kFALSE),
ffdca6aa	106	fCodeDebug(kFALSE),
b65fac7a	107	fJets(0),
04905231	108	fEvent(0),
	109	fTracks(0),
	110	fClus(0)
48d874ff	111	{
48d874ff	112	// Standard constructor.
7efbea04	113
7efbea04	114	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff	115	}
	116
	117	//________________________________________________________________________
914d486c	118	AliEmcalJetTask::~AliEmcalJetTask()
48d874ff	119	{
	120	// Destructor
	121	}
	122
	123	//________________________________________________________________________
914d486c	124	void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff	125	{
	126	// Create user objects.
	127
914d486c	128	fJets = new TClonesArray("AliEmcalJet");
48d874ff	129	fJets->SetName(fJetsName);
48d874ff	130	}
	131
	132	//________________________________________________________________________
914d486c	133	void AliEmcalJetTask::UserExec(Option_t *)
48d874ff	134	{
48d874ff	135	// Main loop, called for each event.
04905231	136	if (!fIsInit) {
	137	if (!DoInit())
	138	return;
	139	fIsInit = kTRUE;
1772fe65	140	}
b65fac7a	141
918dd8c8	142	// clear the jet array (normally a null operation)
	143	fJets->Delete();
	144
04905231	145	FindJets();
48d874ff	146	}
	147
	148	//________________________________________________________________________
914d486c	149	void AliEmcalJetTask::Terminate(Option_t *)
48d874ff	150	{
48d874ff	151	// Called once at the end of the analysis.
48d874ff	152	}
	153
	154	//________________________________________________________________________
04905231	155	void AliEmcalJetTask::FindJets()
48d874ff	156	{
48d874ff	157	// Find jets.
3c5aff5d	158	if (!fTracks && !fClus){
3c5aff5d	159	cout << "WARNING NO TRACKS OR CLUSTERS:" <<endl;
d03084cd	160	return;
3c5aff5d	161	}
3c5aff5d	162
48d874ff	163	TString name("kt");
48d874ff	164	fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
6a20534a	165	if ((fJetType & kAKT) != 0) {
48d874ff	166	name = "antikt";
48d874ff	167	jalgo = fastjet::antikt_algorithm;
ca5c29fa	168	AliDebug(1,"Using AKT algorithm");
	169	}
	170	else {
	171	AliDebug(1,"Using KT algorithm");
48d874ff	172	}
48d874ff	173
6a20534a	174	if ((fJetType & kR020Jet) != 0)
	175	fRadius = 0.2;
	176	else if ((fJetType & kR030Jet) != 0)
	177	fRadius = 0.3;
	178	else if ((fJetType & kR040Jet) != 0)
	179	fRadius = 0.4;
	180
04905231	181	// setup fj wrapper
48d874ff	182	AliFJWrapper fjw(name, name);
f5f3c8e9	183	fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
1f9c287f	184	fjw.SetGhostArea(fGhostArea);
04905231	185	fjw.SetR(fRadius);
7071e730	186	fjw.SetAlgorithm(jalgo);
31bab07e	187	fjw.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
83086e6c	188	fjw.SetMaxRap(fEtaMax);
48d874ff	189	fjw.Clear();
48d874ff	190
04905231	191	// get primary vertex
04905231	192	Double_t vertex[3] = {0, 0, 0};
816bb2c8	193	if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(vertex);
04905231	194
ca5c29fa	195	AliDebug(2,Form("Jet type = %d", fJetType));
6a20534a	196
6a20534a	197	if ((fIsMcPart \|\| ((fJetType & kFullJet) != 0) \|\| ((fJetType & kChargedJet) != 0)) && fTracks) {
04905231	198	const Int_t Ntracks = fTracks->GetEntries();
48d874ff	199	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
04905231	200	AliVParticle t = static_cast<AliVParticle>(fTracks->At(iTracks));
48d874ff	201	if (!t)
48d874ff	202	continue;
6a20534a	203	if (fIsMcPart) {
5be3857d	204	if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
5be3857d	205	AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
6a20534a	206	continue;
5be3857d	207	}
	208	if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
	209	AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
6a20534a	210	continue;
5be3857d	211	}
6a20534a	212	}
a7477843	213	if (fIsMcPart \|\| TMath::Abs(t->GetLabel()) > fMinMCLabel) {
d1f2e0d9	214	if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
d1f2e0d9	215	AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
6a20534a	216	continue;
6a20534a	217	}
d1f2e0d9	218	else {
d1f2e0d9	219	AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
ca5c29fa	220	}
6a20534a	221	}
6a20534a	222	else {
d1f2e0d9	223	if (t->TestBits(fConstSel) != (Int_t)fConstSel) {
d1f2e0d9	224	AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
6a20534a	225	continue;
6a20534a	226	}
d1f2e0d9	227	else {
d1f2e0d9	228	AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
ca5c29fa	229	}
ca5c29fa	230	}
1772fe65	231	if (t->Pt() < fMinJetTrackPt)
914d486c	232	continue;
04905231	233	Double_t eta = t->Eta();
04905231	234	Double_t phi = t->Phi();
ddf74a88	235	if ((eta<fEtaMin) \|\| (eta>fEtaMax) \|\|
ddf74a88	236	(phi<fPhiMin) \|\| (phi>fPhiMax))
04905231	237	continue;
da67b88a	238
6421eeb0	239	// artificial inefficiency
	240	if (fTrackEfficiency < 1.) {
	241	Double_t rnd = gRandom->Rndm();
	242	if (fTrackEfficiency < rnd) {
	243	AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", iTracks));
	244	continue;
	245	}
	246	}
	247
a3c8c8c8	248	// offset of 100 for consistency with cluster ids
ca5c29fa	249	AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt()));
751cb2dc	250	fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
48d874ff	251	}
48d874ff	252	}
6ea168d0	253
6a20534a	254	if ((((fJetType & kFullJet) != 0) \|\| ((fJetType & kNeutralJet) != 0)) && fClus) {
da67b88a	255	const Int_t Nclus = fClus->GetEntries();
	256	for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
	257	AliVCluster *c = 0;
	258	Double_t cEta=0,cPhi=0,cPt=0;
	259	Double_t cPx=0,cPy=0,cPz=0;
	260	if (fIsEmcPart) {
	261	AliEmcalParticle ep = static_cast<AliEmcalParticle>(fClus->At(iClus));
	262	if (!ep)
	263	continue;
6a20534a	264
da67b88a	265	c = ep->GetCluster();
da67b88a	266	if (!c)
6a20534a	267	continue;
6a20534a	268
a7477843	269	if (c->GetLabel() > fMinMCLabel) {
d1f2e0d9	270	if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
d1f2e0d9	271	AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
6a20534a	272	continue;
ca5c29fa	273	}
d1f2e0d9	274	else {
d1f2e0d9	275	AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
ca5c29fa	276	}
6a20534a	277	}
6a20534a	278	else {
d1f2e0d9	279	if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
d1f2e0d9	280	AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
6a20534a	281	continue;
ca5c29fa	282	}
d1f2e0d9	283	else {
d1f2e0d9	284	AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
ca5c29fa	285	}
6a20534a	286	}
6a20534a	287
da67b88a	288	cEta = ep->Eta();
	289	cPhi = ep->Phi();
	290	cPt = ep->Pt();
	291	cPx = ep->Px();
	292	cPy = ep->Py();
	293	cPz = ep->Pz();
	294	} else {
	295	c = static_cast<AliVCluster*>(fClus->At(iClus));
	296	if (!c)
	297	continue;
6a20534a	298
a7477843	299	if (c->GetLabel() > fMinMCLabel) {
d1f2e0d9	300	if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
d1f2e0d9	301	AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
6a20534a	302	continue;
ca5c29fa	303	}
d1f2e0d9	304	else {
d1f2e0d9	305	AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
ca5c29fa	306	}
6a20534a	307	}
6a20534a	308	else {
d1f2e0d9	309	if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
d1f2e0d9	310	AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
6a20534a	311	continue;
ca5c29fa	312	}
d1f2e0d9	313	else {
d1f2e0d9	314	AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
ca5c29fa	315	}
6a20534a	316	}
6a20534a	317
da67b88a	318	TLorentzVector nP;
	319	c->GetMomentum(nP, vertex);
	320	cEta = nP.Eta();
	321	cPhi = nP.Phi();
	322	cPt = nP.Pt();
	323	cPx = nP.Px();
	324	cPy = nP.Py();
	325	cPz = nP.Pz();
04905231	326	}
da67b88a	327	if (!c->IsEMCAL())
	328	continue;
	329	if (cPt < fMinJetClusPt)
	330	continue;
ddf74a88	331	if ((cEta<fEtaMin) \|\| (cEta>fEtaMax) \|\|
ddf74a88	332	(cPhi<fPhiMin) \|\| (cPhi>fPhiMax))
da67b88a	333	continue;
da67b88a	334	// offset of 100 to skip ghost particles uid = -1
ca5c29fa	335	AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
f660c2d6	336	fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPxcPx+cPycPy+cPz*cPz), -iClus - 100);
48d874ff	337	}
48d874ff	338	}
5d0a5007	339
	340	// setting legacy mode
	341	if (fLegacyMode) {
	342	fjw.SetLegacyMode(kTRUE);
	343	}
	344
7efbea04	345	// run jet finder
48d874ff	346	fjw.Run();
f5f3c8e9	347
	348	// get geometry
	349	AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
	350	if (!geom) {
02c7e943	351	AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9	352	return;
f5f3c8e9	353	}
04905231	354
04905231	355	// loop over fastjet jets
48d874ff	356	std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
ca5c29fa	357	// sort jets according to jet pt
	358	static Int_t indexes[9999] = {-1};
	359	GetSortedArray(indexes, jets_incl);
	360
	361	AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
	362	for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
	363	Int_t ij = indexes[ijet];
	364	AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
	365
f5f3c8e9	366	if (jets_incl[ij].perp()<fMinJetPt)
	367	continue;
	368	if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff	369	continue;
c9ad772c	370	if ((jets_incl[ij].eta()<fJetEtaMin) \|\| (jets_incl[ij].eta()>fJetEtaMax) \|\|
	371	(jets_incl[ij].phi()<fJetPhiMin) \|\| (jets_incl[ij].phi()>fJetPhiMax))
	372	continue;
04905231	373
914d486c	374	AliEmcalJet jet = new ((fJets)[jetCount])
914d486c	375	AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
04905231	376
04905231	377	// loop over constituents
c64cb1f6	378	std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
04905231	379	jet->SetNumberOfTracks(constituents.size());
	380	jet->SetNumberOfClusters(constituents.size());
	381
b65fac7a	382	Int_t nt = 0;
b65fac7a	383	Int_t nc = 0;
f5f3c8e9	384	Double_t neutralE = 0;
111318e8	385	Double_t maxCh = 0;
111318e8	386	Double_t maxNe = 0;
f5f3c8e9	387	Int_t gall = 0;
f5f3c8e9	388	Int_t gemc = 0;
04905231	389	Int_t cemc = 0;
111318e8	390	Int_t ncharged = 0;
111318e8	391	Int_t nneutral = 0;
02c7e943	392	Double_t mcpt = 0;
04905231	393	Double_t emcpt = 0;
b65fac7a	394
b65fac7a	395	for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff	396	Int_t uid = constituents[ic].user_index();
5be3857d	397	AliDebug(3,Form("Processing constituent %d", uid));
04905231	398	if ((uid == -1) /&& (constituents[ic].kt2() < 1e-25)/) { //ghost particle
f5f3c8e9	399	++gall;
04905231	400	Double_t gphi = constituents[ic].phi();
	401	if (gphi<0)
	402	gphi += TMath::TwoPi();
	403	gphi *= TMath::RadToDeg();
f5f3c8e9	404	Double_t geta = constituents[ic].eta();
b65fac7a	405	if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
b65fac7a	406	(geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
04905231	407	++gemc;
04905231	408	} else if ((uid > 0) && fTracks) { // track constituent
111318e8	409	Int_t tid = uid - 100;
04905231	410	AliVParticle t = static_cast<AliVParticle>(fTracks->At(tid));
ca5c29fa	411	if (!t) {
ca5c29fa	412	AliError(Form("Could not find track %d",tid));
04905231	413	continue;
ca5c29fa	414	}
	415	if (jetCount < fMarkConst) {
	416	AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
6a20534a	417	t->SetBit(fJetType);
ca5c29fa	418	}
04905231	419	Double_t cEta = t->Eta();
	420	Double_t cPhi = t->Phi();
	421	Double_t cPt = t->Pt();
	422	Double_t cP = t->P();
	423	if (t->Charge() == 0) {
	424	neutralE += cP;
	425	++nneutral;
	426	if (cPt > maxNe)
	427	maxNe = cPt;
	428	} else {
	429	++ncharged;
	430	if (cPt > maxCh)
	431	maxCh = cPt;
	432	}
a7477843	433	if (fIsMcPart \|\| TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
04905231	434	mcpt += cPt;
	435
	436	if (cPhi<0)
	437	cPhi += TMath::TwoPi();
	438	cPhi *= TMath::RadToDeg();
	439	if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
	440	(cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
	441	emcpt += cPt;
	442	++cemc;
d03084cd	443	}
04905231	444
	445	jet->AddTrackAt(tid, nt);
	446	++nt;
	447	} else if (fClus) { // cluster constituent
111318e8	448	Int_t cid = -uid - 100;
04905231	449	AliVCluster *c = 0;
	450	Double_t cEta=0,cPhi=0,cPt=0,cP=0;
	451	if (fIsEmcPart) {
	452	AliEmcalParticle ep = static_cast<AliEmcalParticle>(fClus->At(cid));
	453	if (!ep)
	454	continue;
	455	c = ep->GetCluster();
	456	if (!c)
	457	continue;
7030f36f	458	if (jetCount < fMarkConst)
6a20534a	459	ep->SetBit(fJetType);
04905231	460	cEta = ep->Eta();
	461	cPhi = ep->Phi();
	462	cPt = ep->Pt();
	463	cP = ep->P();
	464	} else {
	465	c = static_cast<AliVCluster*>(fClus->At(cid));
	466	if (!c)
	467	continue;
7030f36f	468	if (jetCount < fMarkConst)
6a20534a	469	c->SetBit(fJetType);
04905231	470	TLorentzVector nP;
	471	c->GetMomentum(nP, vertex);
	472	cEta = nP.Eta();
	473	cPhi = nP.Phi();
	474	cPt = nP.Pt();
	475	cP = nP.P();
	476	}
	477
	478	neutralE += cP;
	479	if (cPt > maxNe)
	480	maxNe = cPt;
	481
a7477843	482	if (c->GetLabel() > fMinMCLabel) // MC particle
43032ce2	483	mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
04905231	484
	485	if (cPhi<0)
	486	cPhi += TMath::TwoPi();
	487	cPhi *= TMath::RadToDeg();
	488	if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
	489	(cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
	490	emcpt += cPt;
	491	++cemc;
	492	}
	493
	494	jet->AddClusterAt(cid, nc);
	495	++nc;
	496	++nneutral;
02c7e943	497	} else {
	498	AliError(Form("%s: No logical way to end up here.", GetName()));
	499	continue;
48d874ff	500	}
04905231	501	} /* end of constituent loop */
04905231	502
35789a2d	503	jet->SetNumberOfTracks(nt);
35789a2d	504	jet->SetNumberOfClusters(nc);
f5f3c8e9	505	jet->SortConstituents();
111318e8	506	jet->SetMaxNeutralPt(maxNe);
111318e8	507	jet->SetMaxChargedPt(maxCh);
b65fac7a	508	jet->SetNEF(neutralE / jet->E());
a88b5c99	509	fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
	510	jet->SetArea(area.perp());
	511	jet->SetAreaEta(area.eta());
	512	jet->SetAreaPhi(area.phi());
111318e8	513	jet->SetNumberOfCharged(ncharged);
111318e8	514	jet->SetNumberOfNeutrals(nneutral);
02c7e943	515	jet->SetMCPt(mcpt);
04905231	516	jet->SetNEmc(cemc);
	517	jet->SetPtEmc(emcpt);
	518
b65fac7a	519	if (gall > 0)
b65fac7a	520	jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9	521	else
f5f3c8e9	522	jet->SetAreaEmc(-1);
b65fac7a	523	if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65	524	(jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
	525	(jet->Eta() > geom->GetArm1EtaMin()) &&
	526	(jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9	527	jet->SetAxisInEmcal(kTRUE);
ca5c29fa	528
ca5c29fa	529	AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
48d874ff	530	jetCount++;
48d874ff	531	}
7030f36f	532	//fJets->Sort();
48d874ff	533	}
04905231	534
ca5c29fa	535	//________________________________________________________________________
	536	Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
	537	{
	538	// Get the leading jets.
	539
	540	static Float_t pt[9999] = {0};
	541
	542	const Int_t n = (Int_t)array.size();
	543
	544	if (n < 1)
	545	return kFALSE;
	546
	547	for (Int_t i = 0; i < n; i++)
	548	pt[i] = array[i].perp();
	549
	550	TMath::Sort(n, pt, indexes);
	551
	552	return kTRUE;
	553	}
	554
04905231	555	//________________________________________________________________________
	556	Bool_t AliEmcalJetTask::DoInit()
	557	{
	558	// Init. Return true if successful.
	559
6421eeb0	560	if (fTrackEfficiency < 1.) {
	561	if (gRandom) delete gRandom;
	562	gRandom = new TRandom3(0);
	563	}
	564
04905231	565	// get input collections
	566	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	567
	568	// get the event
	569	fEvent = InputEvent();
	570	if (!fEvent) {
	571	AliError(Form("%s: Could not retrieve event! Returning", GetName()));
	572	return 0;
	573	}
	574
	575	// add jets to event if not yet there
	576	if (!(fEvent->FindListObject(fJetsName)))
	577	fEvent->AddObject(fJets);
	578	else {
f4b49fa6	579	AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
04905231	580	return 0;
	581	}
	582
da67b88a	583	if (fTracksName == "Tracks")
	584	am->LoadBranch("Tracks");
	585	if (!fTracks && !fTracksName.IsNull()) {
	586	fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
	587	if (!fTracks) {
	588	AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
	589	return 0;
04905231	590	}
da67b88a	591	}
	592	if (fTracks) {
	593	TClass cls(fTracks->GetClass()->GetName());
5be3857d	594	if (cls.InheritsFrom("AliMCParticle") \|\| cls.InheritsFrom("AliAODMCParticle"))
04905231	595	fIsMcPart = 1;
04905231	596	}
da67b88a	597
	598	if (fCaloName == "CaloClusters")
	599	am->LoadBranch("CaloClusters");
	600	if (!fClus && !fCaloName.IsNull()) {
	601	fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
	602	if (!fClus) {
	603	AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
	604	return 0;
04905231	605	}
da67b88a	606	}
	607	if (fClus) {
	608	TClass cls(fClus->GetClass()->GetName());
04905231	609	if (cls.InheritsFrom("AliEmcalParticle"))
	610	fIsEmcPart = 1;
	611	}
	612
	613	return 1;
	614	}