[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::BIpt_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::BIpt_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()));
      if (fCodeDebug) { 
        fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
      } else {
        fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), 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),
31bab07e	59	fRecombScheme(fastjet::BIpt_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),
31bab07e	100	fRecombScheme(fastjet::BIpt_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()));
31bab07e	251	if (fCodeDebug) {
ffdca6aa	252	fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
31bab07e	253	} else {
ffdca6aa	254	fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
31bab07e	255	}
48d874ff	256	}
48d874ff	257	}
6ea168d0	258
6a20534a	259	if ((((fJetType & kFullJet) != 0) \|\| ((fJetType & kNeutralJet) != 0)) && fClus) {
da67b88a	260	const Int_t Nclus = fClus->GetEntries();
	261	for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
	262	AliVCluster *c = 0;
	263	Double_t cEta=0,cPhi=0,cPt=0;
	264	Double_t cPx=0,cPy=0,cPz=0;
	265	if (fIsEmcPart) {
	266	AliEmcalParticle ep = static_cast<AliEmcalParticle>(fClus->At(iClus));
	267	if (!ep)
	268	continue;
6a20534a	269
da67b88a	270	c = ep->GetCluster();
da67b88a	271	if (!c)
6a20534a	272	continue;
6a20534a	273
a7477843	274	if (c->GetLabel() > fMinMCLabel) {
d1f2e0d9	275	if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
d1f2e0d9	276	AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
6a20534a	277	continue;
ca5c29fa	278	}
d1f2e0d9	279	else {
d1f2e0d9	280	AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
ca5c29fa	281	}
6a20534a	282	}
6a20534a	283	else {
d1f2e0d9	284	if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
d1f2e0d9	285	AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
6a20534a	286	continue;
ca5c29fa	287	}
d1f2e0d9	288	else {
d1f2e0d9	289	AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
ca5c29fa	290	}
6a20534a	291	}
6a20534a	292
da67b88a	293	cEta = ep->Eta();
	294	cPhi = ep->Phi();
	295	cPt = ep->Pt();
	296	cPx = ep->Px();
	297	cPy = ep->Py();
	298	cPz = ep->Pz();
	299	} else {
	300	c = static_cast<AliVCluster*>(fClus->At(iClus));
	301	if (!c)
	302	continue;
6a20534a	303
a7477843	304	if (c->GetLabel() > fMinMCLabel) {
d1f2e0d9	305	if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
d1f2e0d9	306	AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
6a20534a	307	continue;
ca5c29fa	308	}
d1f2e0d9	309	else {
d1f2e0d9	310	AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
ca5c29fa	311	}
6a20534a	312	}
6a20534a	313	else {
d1f2e0d9	314	if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
d1f2e0d9	315	AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
6a20534a	316	continue;
ca5c29fa	317	}
d1f2e0d9	318	else {
d1f2e0d9	319	AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
ca5c29fa	320	}
6a20534a	321	}
6a20534a	322
da67b88a	323	TLorentzVector nP;
	324	c->GetMomentum(nP, vertex);
	325	cEta = nP.Eta();
	326	cPhi = nP.Phi();
	327	cPt = nP.Pt();
	328	cPx = nP.Px();
	329	cPy = nP.Py();
	330	cPz = nP.Pz();
04905231	331	}
da67b88a	332	if (!c->IsEMCAL())
	333	continue;
	334	if (cPt < fMinJetClusPt)
	335	continue;
ddf74a88	336	if ((cEta<fEtaMin) \|\| (cEta>fEtaMax) \|\|
ddf74a88	337	(cPhi<fPhiMin) \|\| (cPhi>fPhiMax))
da67b88a	338	continue;
da67b88a	339	// offset of 100 to skip ghost particles uid = -1
ca5c29fa	340	AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
f660c2d6	341	fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPxcPx+cPycPy+cPz*cPz), -iClus - 100);
48d874ff	342	}
48d874ff	343	}
5d0a5007	344
	345	// setting legacy mode
	346	if (fLegacyMode) {
	347	fjw.SetLegacyMode(kTRUE);
	348	}
	349
7efbea04	350	// run jet finder
48d874ff	351	fjw.Run();
f5f3c8e9	352
	353	// get geometry
	354	AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
	355	if (!geom) {
02c7e943	356	AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9	357	return;
f5f3c8e9	358	}
04905231	359
04905231	360	// loop over fastjet jets
48d874ff	361	std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
ca5c29fa	362	// sort jets according to jet pt
	363	static Int_t indexes[9999] = {-1};
	364	GetSortedArray(indexes, jets_incl);
	365
	366	AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
	367	for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
	368	Int_t ij = indexes[ijet];
	369	AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
	370
f5f3c8e9	371	if (jets_incl[ij].perp()<fMinJetPt)
	372	continue;
	373	if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff	374	continue;
c9ad772c	375	if ((jets_incl[ij].eta()<fJetEtaMin) \|\| (jets_incl[ij].eta()>fJetEtaMax) \|\|
	376	(jets_incl[ij].phi()<fJetPhiMin) \|\| (jets_incl[ij].phi()>fJetPhiMax))
	377	continue;
04905231	378
914d486c	379	AliEmcalJet jet = new ((fJets)[jetCount])
914d486c	380	AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
04905231	381
04905231	382	// loop over constituents
c64cb1f6	383	std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
04905231	384	jet->SetNumberOfTracks(constituents.size());
	385	jet->SetNumberOfClusters(constituents.size());
	386
b65fac7a	387	Int_t nt = 0;
b65fac7a	388	Int_t nc = 0;
f5f3c8e9	389	Double_t neutralE = 0;
111318e8	390	Double_t maxCh = 0;
111318e8	391	Double_t maxNe = 0;
f5f3c8e9	392	Int_t gall = 0;
f5f3c8e9	393	Int_t gemc = 0;
04905231	394	Int_t cemc = 0;
111318e8	395	Int_t ncharged = 0;
111318e8	396	Int_t nneutral = 0;
02c7e943	397	Double_t mcpt = 0;
04905231	398	Double_t emcpt = 0;
b65fac7a	399
b65fac7a	400	for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff	401	Int_t uid = constituents[ic].user_index();
5be3857d	402	AliDebug(3,Form("Processing constituent %d", uid));
04905231	403	if ((uid == -1) /&& (constituents[ic].kt2() < 1e-25)/) { //ghost particle
f5f3c8e9	404	++gall;
04905231	405	Double_t gphi = constituents[ic].phi();
	406	if (gphi<0)
	407	gphi += TMath::TwoPi();
	408	gphi *= TMath::RadToDeg();
f5f3c8e9	409	Double_t geta = constituents[ic].eta();
b65fac7a	410	if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
b65fac7a	411	(geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
04905231	412	++gemc;
04905231	413	} else if ((uid > 0) && fTracks) { // track constituent
111318e8	414	Int_t tid = uid - 100;
04905231	415	AliVParticle t = static_cast<AliVParticle>(fTracks->At(tid));
ca5c29fa	416	if (!t) {
ca5c29fa	417	AliError(Form("Could not find track %d",tid));
04905231	418	continue;
ca5c29fa	419	}
	420	if (jetCount < fMarkConst) {
	421	AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
6a20534a	422	t->SetBit(fJetType);
ca5c29fa	423	}
04905231	424	Double_t cEta = t->Eta();
	425	Double_t cPhi = t->Phi();
	426	Double_t cPt = t->Pt();
	427	Double_t cP = t->P();
	428	if (t->Charge() == 0) {
	429	neutralE += cP;
	430	++nneutral;
	431	if (cPt > maxNe)
	432	maxNe = cPt;
	433	} else {
	434	++ncharged;
	435	if (cPt > maxCh)
	436	maxCh = cPt;
	437	}
a7477843	438	if (fIsMcPart \|\| TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
04905231	439	mcpt += cPt;
	440
	441	if (cPhi<0)
	442	cPhi += TMath::TwoPi();
	443	cPhi *= TMath::RadToDeg();
	444	if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
	445	(cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
	446	emcpt += cPt;
	447	++cemc;
d03084cd	448	}
04905231	449
	450	jet->AddTrackAt(tid, nt);
	451	++nt;
	452	} else if (fClus) { // cluster constituent
111318e8	453	Int_t cid = -uid - 100;
04905231	454	AliVCluster *c = 0;
	455	Double_t cEta=0,cPhi=0,cPt=0,cP=0;
	456	if (fIsEmcPart) {
	457	AliEmcalParticle ep = static_cast<AliEmcalParticle>(fClus->At(cid));
	458	if (!ep)
	459	continue;
	460	c = ep->GetCluster();
	461	if (!c)
	462	continue;
7030f36f	463	if (jetCount < fMarkConst)
6a20534a	464	ep->SetBit(fJetType);
04905231	465	cEta = ep->Eta();
	466	cPhi = ep->Phi();
	467	cPt = ep->Pt();
	468	cP = ep->P();
	469	} else {
	470	c = static_cast<AliVCluster*>(fClus->At(cid));
	471	if (!c)
	472	continue;
7030f36f	473	if (jetCount < fMarkConst)
6a20534a	474	c->SetBit(fJetType);
04905231	475	TLorentzVector nP;
	476	c->GetMomentum(nP, vertex);
	477	cEta = nP.Eta();
	478	cPhi = nP.Phi();
	479	cPt = nP.Pt();
	480	cP = nP.P();
	481	}
	482
	483	neutralE += cP;
	484	if (cPt > maxNe)
	485	maxNe = cPt;
	486
a7477843	487	if (c->GetLabel() > fMinMCLabel) // MC particle
43032ce2	488	mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
04905231	489
	490	if (cPhi<0)
	491	cPhi += TMath::TwoPi();
	492	cPhi *= TMath::RadToDeg();
	493	if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
	494	(cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
	495	emcpt += cPt;
	496	++cemc;
	497	}
	498
	499	jet->AddClusterAt(cid, nc);
	500	++nc;
	501	++nneutral;
02c7e943	502	} else {
	503	AliError(Form("%s: No logical way to end up here.", GetName()));
	504	continue;
48d874ff	505	}
04905231	506	} /* end of constituent loop */
04905231	507
35789a2d	508	jet->SetNumberOfTracks(nt);
35789a2d	509	jet->SetNumberOfClusters(nc);
f5f3c8e9	510	jet->SortConstituents();
111318e8	511	jet->SetMaxNeutralPt(maxNe);
111318e8	512	jet->SetMaxChargedPt(maxCh);
b65fac7a	513	jet->SetNEF(neutralE / jet->E());
a88b5c99	514	fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
	515	jet->SetArea(area.perp());
	516	jet->SetAreaEta(area.eta());
	517	jet->SetAreaPhi(area.phi());
111318e8	518	jet->SetNumberOfCharged(ncharged);
111318e8	519	jet->SetNumberOfNeutrals(nneutral);
02c7e943	520	jet->SetMCPt(mcpt);
04905231	521	jet->SetNEmc(cemc);
	522	jet->SetPtEmc(emcpt);
	523
b65fac7a	524	if (gall > 0)
b65fac7a	525	jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9	526	else
f5f3c8e9	527	jet->SetAreaEmc(-1);
b65fac7a	528	if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65	529	(jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
	530	(jet->Eta() > geom->GetArm1EtaMin()) &&
	531	(jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9	532	jet->SetAxisInEmcal(kTRUE);
ca5c29fa	533
ca5c29fa	534	AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
48d874ff	535	jetCount++;
48d874ff	536	}
7030f36f	537	//fJets->Sort();
48d874ff	538	}
04905231	539
ca5c29fa	540	//________________________________________________________________________
	541	Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
	542	{
	543	// Get the leading jets.
	544
	545	static Float_t pt[9999] = {0};
	546
	547	const Int_t n = (Int_t)array.size();
	548
	549	if (n < 1)
	550	return kFALSE;
	551
	552	for (Int_t i = 0; i < n; i++)
	553	pt[i] = array[i].perp();
	554
	555	TMath::Sort(n, pt, indexes);
	556
	557	return kTRUE;
	558	}
	559
04905231	560	//________________________________________________________________________
	561	Bool_t AliEmcalJetTask::DoInit()
	562	{
	563	// Init. Return true if successful.
	564
6421eeb0	565	if (fTrackEfficiency < 1.) {
	566	if (gRandom) delete gRandom;
	567	gRandom = new TRandom3(0);
	568	}
	569
04905231	570	// get input collections
	571	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	572
	573	// get the event
	574	fEvent = InputEvent();
	575	if (!fEvent) {
	576	AliError(Form("%s: Could not retrieve event! Returning", GetName()));
	577	return 0;
	578	}
	579
	580	// add jets to event if not yet there
	581	if (!(fEvent->FindListObject(fJetsName)))
	582	fEvent->AddObject(fJets);
	583	else {
f4b49fa6	584	AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
04905231	585	return 0;
	586	}
	587
da67b88a	588	if (fTracksName == "Tracks")
	589	am->LoadBranch("Tracks");
	590	if (!fTracks && !fTracksName.IsNull()) {
	591	fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
	592	if (!fTracks) {
	593	AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
	594	return 0;
04905231	595	}
da67b88a	596	}
	597	if (fTracks) {
	598	TClass cls(fTracks->GetClass()->GetName());
5be3857d	599	if (cls.InheritsFrom("AliMCParticle") \|\| cls.InheritsFrom("AliAODMCParticle"))
04905231	600	fIsMcPart = 1;
04905231	601	}
da67b88a	602
	603	if (fCaloName == "CaloClusters")
	604	am->LoadBranch("CaloClusters");
	605	if (!fClus && !fCaloName.IsNull()) {
	606	fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
	607	if (!fClus) {
	608	AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
	609	return 0;
04905231	610	}
da67b88a	611	}
	612	if (fClus) {
	613	TClass cls(fClus->GetClass()->GetName());
04905231	614	if (cls.InheritsFrom("AliEmcalParticle"))
	615	fIsEmcPart = 1;
	616	}
	617
	618	return 1;
	619	}