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

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