[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),
  fTrackEfficiency(1.),
  fIsInit(0),
  fIsPSelSet(0),
  fIsMcPart(0),
  fIsEmcPart(0),
  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),
  fTrackEfficiency(1.),
  fIsInit(0),
  fIsPSelSet(0),
  fIsMcPart(0),
  fIsEmcPart(0),
  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;
  }

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