[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"

ClassImp(AliEmcalJetTask)

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask() : 
  AliAnalysisTaskSE("AliEmcalJetTask"),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fJetType(kNone),
  fConstSel(kAllJets),
  fMCConstSel(kAllJets),
  fMarkConst(kFALSE),
  fRadius(0.4),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fPhiMin(0),
  fPhiMax(TMath::TwoPi()),
  fEtaMin(-1),
  fEtaMax(1),
  fMinJetArea(0.01),
  fMinJetPt(1.0),
  fGhostArea(0.01),
  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(kAllJets),
  fMCConstSel(kAllJets),
  fMarkConst(kFALSE),
  fRadius(0.4),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fPhiMin(0),
  fPhiMax(TMath::TwoPi()),
  fEtaMin(-1),
  fEtaMax(1),
  fMinJetArea(0.01),
  fMinJetPt(1.0),
  fGhostArea(0.01),
  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)
    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(1);
  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))
	  continue;
	if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0))
	  continue;
      }
      if (fIsMcPart || t->GetLabel() != 0) {
	if (fMCConstSel == kNone) {
	  AliDebug(2,Form("Skipping track %d because bit mask is 0.", iTracks));
	  continue;
	}
	if (fMCConstSel != kAllJets) {
	  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 (fConstSel == kNone) {
	  AliDebug(2,Form("Skipping track %d because bit mask is 0.", iTracks));
	  continue;
	}
	if (fConstSel != kAllJets) {
	  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() > 0) {
	  if (fMCConstSel == kNone) {
	    AliDebug(2,Form("Skipping cluster %d because bit mask is 0.", iClus));
	    continue;
	  }
	  if (fMCConstSel != kAllJets) {
	    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 (fConstSel == kNone) {
	    AliDebug(2,Form("Skipping cluster %d because bit mask is 0.", iClus));
	    continue;
	  }
	  if (fConstSel != kAllJets) {
	    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() > 0) {
	  if (fMCConstSel == kNone) {
	    AliDebug(2,Form("Skipping cluster %d because bit mask is 0.", iClus));
	    continue;
	  }
	  if (fMCConstSel != kAllJets) {
	    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 (fConstSel == kNone) {
	    AliDebug(2,Form("Skipping cluster %d because bit mask is 0.", iClus));
	    continue;
	  }
	  if (fConstSel != kAllJets) {
	    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;

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