[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 "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),
  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),
  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;

      // 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 += cPt * c->GetMCEnergyFraction();

        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.

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