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

// $Id$
//
// Emcal jet finder task.
//
// Authors: C.Loizides, S.Aiola

#include <TChain.h>
#include <TClonesArray.h>
#include <TList.h>
#include <TLorentzVector.h>

#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliEMCALGeometry.h"
#include "AliEmcalJet.h"
#include "AliFJWrapper.h"
#include "AliVCluster.h"
#include "AliVParticle.h"
#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"

#include "AliEmcalJetTask.h"

ClassImp(AliEmcalJetTask)

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask() : 
  AliAnalysisTaskSE("AliEmcalJetTask"),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fAlgo(1),
  fRadius(0.4),
  fType(0),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fMinJetArea(0.01),
  fMinJetPt(1.0),
  fJets(0),
  fEvent(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask(const char *name) : 
  AliAnalysisTaskSE(name),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fAlgo(1),
  fRadius(0.4),
  fType(0),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fMinJetArea(0.01),
  fMinJetPt(1.0),
  fJets(0),
  fEvent(0)
{
  // Standard constructor.

  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
}

//________________________________________________________________________
AliEmcalJetTask::~AliEmcalJetTask()
{
  // Destructor
}

//________________________________________________________________________
void AliEmcalJetTask::UserCreateOutputObjects()
{
  // Create user objects.

  fJets = new TClonesArray("AliEmcalJet");
  fJets->SetName(fJetsName);
}

//_____________________________________________________
TString AliEmcalJetTask::GetBeamType()
{
  // Get beam type : pp-AA-pA
  // ESDs have it directly, AODs get it from hardcoded run number ranges

  TString beamType;

  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fEvent);
  if (esd) {
    const AliESDRun *run = esd->GetESDRun();
    beamType = run->GetBeamType();
  }
  else
  {
    Int_t runNumber = fEvent->GetRunNumber();
    if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
	(runNumber >= 166529 && runNumber <= 170593))    // LHC11h
    {
      beamType = "A-A";
    }
    else 
    {
      beamType = "p-p";
    }
  }

  return beamType;    
}


//________________________________________________________________________
void AliEmcalJetTask::UserExec(Option_t *) 
{
  // Main loop, called for each event.

  // get the event
  fEvent = InputEvent();

  if (!fEvent) {
    AliError(Form("%s: Could not retrieve event! Returning", GetName()));
    return;
  }

  // add jets to event if not yet there
  if (!(fEvent->FindListObject(fJetsName)))
    fEvent->AddObject(fJets);

  // delete jet output
  fJets->Delete();

  // get input collections
  TClonesArray *tracks = 0;
  TClonesArray *clus   = 0;
  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();

  if ((fType == 0) || (fType == 1)) {
    if (fTracksName == "Tracks")
      am->LoadBranch("Tracks");
    tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
    if (!tracks) {
      AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
      return;
    }
  }
  if ((fType == 0) || (fType == 2)) {
    if (fCaloName == "CaloClusters")
      am->LoadBranch("CaloClusters");
    clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
    if (!clus) {
      AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
      return;
    }
  }
  
  // get centrality 
  Double_t cent = 99; 
  
  if (GetBeamType() == "A-A") {
    AliCentrality *centrality = InputEvent()->GetCentrality();
    
    if (centrality)
      cent = centrality->GetCentralityPercentile("V0M");
    else
      cent = 99; // probably pp data
    
    if (cent < 0) {
      AliWarning(Form("Centrality negative: %f, assuming 99", cent));
      cent = 99;
    }
  }

  FindJets(tracks, clus, fAlgo, fRadius, cent);
}

//________________________________________________________________________
void AliEmcalJetTask::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}

//________________________________________________________________________
void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
{
  // Find jets.

  if (!tracks && !clus)
    return;

  TString name("kt");
  fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
  if (algo>=1) {
    name  = "antikt";
    jalgo = fastjet::antikt_algorithm;
  }

  AliFJWrapper fjw(name, name);
  fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
  fjw.SetR(radius);
  fjw.SetAlgorithm(jalgo);  
  fjw.SetMaxRap(0.9);
  fjw.Clear();

  if (tracks) {
    const Int_t Ntracks = tracks->GetEntries();
    for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
      AliVParticle *t = dynamic_cast<AliVParticle*>(tracks->At(iTracks));
      if (!t)
        continue;

      if (t->Pt() < fMinJetTrackPt) 
        continue;
      // offset of 100 for consistency with cluster ids
      fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);  
    }
  }

  if (clus) {
    Double_t vertex[3] = {0, 0, 0};
    fEvent->GetPrimaryVertex()->GetXYZ(vertex);
    const Int_t Nclus = clus->GetEntries();
    for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
      AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
      if (!c)
        continue;
      if (!c->IsEMCAL())
        continue;
      TLorentzVector nPart;
      c->GetMomentum(nPart, vertex);
      if (nPart.Pt() < fMinJetClusPt) 
        continue;
      // offset of 100 to skip ghost particles uid = -1
      fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);  
    }
  }

  // run jet finder
  fjw.Run();

  // get geometry
  AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
  if (!geom) {
    AliFatal(Form("%s: Can not create geometry", GetName()));
    return;
  }
  
  std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
  for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++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());
    Double_t vertex[3] = {0, 0, 0};
    fEvent->GetPrimaryVertex()->GetXYZ(vertex);
    vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
    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 ncharged      = 0;
    Int_t nneutral      = 0;
    Double_t mcpt       = 0;

    jet->SetNumberOfTracks(constituents.size());
    jet->SetNumberOfClusters(constituents.size());
    for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
      Int_t uid = constituents[ic].user_index();

      if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
        ++gall;
        Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
        Double_t geta = constituents[ic].eta();
        if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
            (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
          ++gemc;
        continue;
      }	else if ((uid > 0) && tracks) { // track constituent
	Int_t tid = uid - 100;
        AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
        if (t) {
	  if (t->Charge() == 0) {
	    neutralE += t->P();
	    ++nneutral;
	    if (t->Pt() > maxNe)
	      maxNe = t->Pt();
	  } else {
	    ++ncharged;
	    if (t->Pt() > maxCh)
	      maxCh = t->Pt();
	  }

	  if (t->InheritsFrom("AliMCParticle") || t->GetLabel() == 100) // MC particle
	    mcpt += t->Pt();

          jet->AddTrackAt(tid, nt);
          ++nt;
        }
      } else if (clus) { // cluster constituent
	Int_t cid = -uid - 100;
	AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
	if (c) {
	  TLorentzVector nP;
	  c->GetMomentum(nP, vertex);
	  neutralE += nP.P();
	  if (nP.Pt() > maxNe)
	    maxNe = nP.Pt();

	  if (c->Chi2() == 100) // MC particle
	    mcpt += nP.Pt();

	  jet->AddClusterAt(cid, nc);
	  ++nc;
	  ++nneutral;
	}
      } else {
        AliError(Form("%s: No logical way to end up here.", GetName()));
        continue;
      }
    }
    jet->SetNumberOfTracks(nt);
    jet->SetNumberOfClusters(nc);
    jet->SortConstituents();
    jet->SetMaxNeutralPt(maxNe);
    jet->SetMaxChargedPt(maxCh);
    jet->SetNEF(neutralE / jet->E());
    jet->SetArea(fjw.GetJetArea(ij));
    jet->SetNumberOfCharged(ncharged);
    jet->SetNumberOfNeutrals(nneutral);
    jet->SetMCPt(mcpt);
    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);
    jetCount++;
  }
}
Commit	Line	Data
48d874ff	1	// $Id$
542ac102	2	//
	3	// Emcal jet finder task.
	4	//
b65fac7a	5	// Authors: C.Loizides, S.Aiola
48d874ff	6
6ea168d0	7	#include <TChain.h>
48d874ff	8	#include <TClonesArray.h>
6ea168d0	9	#include <TList.h>
6ea168d0	10	#include <TLorentzVector.h>
b65fac7a	11
48d874ff	12	#include "AliAnalysisManager.h"
6ea168d0	13	#include "AliCentrality.h"
f5f3c8e9	14	#include "AliEMCALGeometry.h"
f5f3c8e9	15	#include "AliEmcalJet.h"
63206072	16	#include "AliFJWrapper.h"
35789a2d	17	#include "AliVCluster.h"
b65fac7a	18	#include "AliVParticle.h"
b65fac7a	19	#include "AliVEvent.h"
111318e8	20	#include "AliESDEvent.h"
b65fac7a	21	#include "AliMCEvent.h"
	22
	23	#include "AliEmcalJetTask.h"
48d874ff	24
914d486c	25	ClassImp(AliEmcalJetTask)
48d874ff	26
6ea168d0	27	//________________________________________________________________________
914d486c	28	AliEmcalJetTask::AliEmcalJetTask() :
914d486c	29	AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0	30	fTracksName("Tracks"),
	31	fCaloName("CaloClusters"),
	32	fJetsName("Jets"),
	33	fAlgo(1),
	34	fRadius(0.4),
	35	fType(0),
914d486c	36	fMinJetTrackPt(0.15),
914d486c	37	fMinJetClusPt(0.15),
f5f3c8e9	38	fMinJetArea(0.01),
f5f3c8e9	39	fMinJetPt(1.0),
b65fac7a	40	fJets(0),
b65fac7a	41	fEvent(0)
6ea168d0	42	{
6ea168d0	43	// Default constructor.
6ea168d0	44	}
6ea168d0	45
48d874ff	46	//________________________________________________________________________
914d486c	47	AliEmcalJetTask::AliEmcalJetTask(const char *name) :
63206072	48	AliAnalysisTaskSE(name),
7efbea04	49	fTracksName("Tracks"),
48d874ff	50	fCaloName("CaloClusters"),
7efbea04	51	fJetsName("Jets"),
48d874ff	52	fAlgo(1),
	53	fRadius(0.4),
	54	fType(0),
6ea168d0	55	fMinJetTrackPt(0.15),
914d486c	56	fMinJetClusPt(0.15),
f5f3c8e9	57	fMinJetArea(0.01),
f5f3c8e9	58	fMinJetPt(1.0),
b65fac7a	59	fJets(0),
b65fac7a	60	fEvent(0)
48d874ff	61	{
48d874ff	62	// Standard constructor.
7efbea04	63
7efbea04	64	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff	65	}
	66
	67	//________________________________________________________________________
914d486c	68	AliEmcalJetTask::~AliEmcalJetTask()
48d874ff	69	{
	70	// Destructor
	71	}
	72
	73	//________________________________________________________________________
914d486c	74	void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff	75	{
	76	// Create user objects.
	77
914d486c	78	fJets = new TClonesArray("AliEmcalJet");
48d874ff	79	fJets->SetName(fJetsName);
48d874ff	80	}
48d874ff	81
111318e8	82	//_____________________________________________________
	83	TString AliEmcalJetTask::GetBeamType()
	84	{
	85	// Get beam type : pp-AA-pA
	86	// ESDs have it directly, AODs get it from hardcoded run number ranges
	87
	88	TString beamType;
	89
	90	AliESDEvent esd = dynamic_cast<AliESDEvent>(fEvent);
	91	if (esd) {
	92	const AliESDRun *run = esd->GetESDRun();
	93	beamType = run->GetBeamType();
	94	}
	95	else
	96	{
	97	Int_t runNumber = fEvent->GetRunNumber();
	98	if ((runNumber >= 136851 && runNumber <= 139517) \|\| // LHC10h
	99	(runNumber >= 166529 && runNumber <= 170593)) // LHC11h
	100	{
	101	beamType = "A-A";
	102	}
	103	else
	104	{
	105	beamType = "p-p";
	106	}
	107	}
	108
	109	return beamType;
	110	}
	111
	112
48d874ff	113	//________________________________________________________________________
914d486c	114	void AliEmcalJetTask::UserExec(Option_t *)
48d874ff	115	{
	116	// Main loop, called for each event.
	117
b65fac7a	118	// get the event
1772fe65	119	fEvent = InputEvent();
	120
	121	if (!fEvent) {
02c7e943	122	AliError(Form("%s: Could not retrieve event! Returning", GetName()));
1772fe65	123	return;
1772fe65	124	}
b65fac7a	125
63206072	126	// add jets to event if not yet there
b65fac7a	127	if (!(fEvent->FindListObject(fJetsName)))
b65fac7a	128	fEvent->AddObject(fJets);
48d874ff	129
63206072	130	// delete jet output
	131	fJets->Delete();
	132
63206072	133	// get input collections
48d874ff	134	TClonesArray *tracks = 0;
48d874ff	135	TClonesArray *clus = 0;
63206072	136	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
6ea168d0	137
b65fac7a	138	if ((fType == 0) \|\| (fType == 1)) {
7efbea04	139	if (fTracksName == "Tracks")
48d874ff	140	am->LoadBranch("Tracks");
b65fac7a	141	tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
48d874ff	142	if (!tracks) {
02c7e943	143	AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
48d874ff	144	return;
	145	}
	146	}
b65fac7a	147	if ((fType == 0) \|\| (fType == 2)) {
48d874ff	148	if (fCaloName == "CaloClusters")
48d874ff	149	am->LoadBranch("CaloClusters");
b65fac7a	150	clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
48d874ff	151	if (!clus) {
02c7e943	152	AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
48d874ff	153	return;
48d874ff	154	}
48d874ff	155	}
111318e8	156
	157	// get centrality
	158	Double_t cent = 99;
	159
	160	if (GetBeamType() == "A-A") {
	161	AliCentrality *centrality = InputEvent()->GetCentrality();
	162
	163	if (centrality)
	164	cent = centrality->GetCentralityPercentile("V0M");
	165	else
	166	cent = 99; // probably pp data
	167
	168	if (cent < 0) {
	169	AliWarning(Form("Centrality negative: %f, assuming 99", cent));
	170	cent = 99;
	171	}
63206072	172	}
653159b9	173
6ea168d0	174	FindJets(tracks, clus, fAlgo, fRadius, cent);
48d874ff	175	}
	176
	177	//________________________________________________________________________
914d486c	178	void AliEmcalJetTask::Terminate(Option_t *)
48d874ff	179	{
48d874ff	180	// Called once at the end of the analysis.
48d874ff	181	}
	182
	183	//________________________________________________________________________
914d486c	184	void AliEmcalJetTask::FindJets(TObjArray tracks, TObjArray clus, Int_t algo, Double_t radius, Float_t /cent/)
48d874ff	185	{
	186	// Find jets.
	187
d03084cd	188	if (!tracks && !clus)
	189	return;
	190
48d874ff	191	TString name("kt");
	192	fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
	193	if (algo>=1) {
	194	name = "antikt";
	195	jalgo = fastjet::antikt_algorithm;
	196	}
	197
	198	AliFJWrapper fjw(name, name);
f5f3c8e9	199	fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
48d874ff	200	fjw.SetR(radius);
6ea168d0	201	fjw.SetAlgorithm(jalgo);
48d874ff	202	fjw.SetMaxRap(0.9);
	203	fjw.Clear();
	204
	205	if (tracks) {
	206	const Int_t Ntracks = tracks->GetEntries();
	207	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1772fe65	208	AliVParticle t = dynamic_cast<AliVParticle>(tracks->At(iTracks));
48d874ff	209	if (!t)
48d874ff	210	continue;
1772fe65	211
1772fe65	212	if (t->Pt() < fMinJetTrackPt)
914d486c	213	continue;
a3c8c8c8	214	// offset of 100 for consistency with cluster ids
a3c8c8c8	215	fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
48d874ff	216	}
48d874ff	217	}
6ea168d0	218
48d874ff	219	if (clus) {
48d874ff	220	Double_t vertex[3] = {0, 0, 0};
b65fac7a	221	fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff	222	const Int_t Nclus = clus->GetEntries();
914d486c	223	for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
7efbea04	224	AliVCluster c = dynamic_cast<AliVCluster>(clus->At(iClus));
d03084cd	225	if (!c)
d03084cd	226	continue;
48d874ff	227	if (!c->IsEMCAL())
	228	continue;
	229	TLorentzVector nPart;
	230	c->GetMomentum(nPart, vertex);
111318e8	231	if (nPart.Pt() < fMinJetClusPt)
914d486c	232	continue;
a3c8c8c8	233	// offset of 100 to skip ghost particles uid = -1
a3c8c8c8	234	fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);
48d874ff	235	}
	236	}
	237
7efbea04	238	// run jet finder
48d874ff	239	fjw.Run();
f5f3c8e9	240
	241	// get geometry
	242	AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
	243	if (!geom) {
02c7e943	244	AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9	245	return;
f5f3c8e9	246	}
6ea168d0	247
48d874ff	248	std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
b65fac7a	249	for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9	250	if (jets_incl[ij].perp()<fMinJetPt)
	251	continue;
	252	if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff	253	continue;
914d486c	254	AliEmcalJet jet = new ((fJets)[jetCount])
914d486c	255	AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
6ea168d0	256	Double_t vertex[3] = {0, 0, 0};
b65fac7a	257	fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff	258	vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
b65fac7a	259	Int_t nt = 0;
b65fac7a	260	Int_t nc = 0;
f5f3c8e9	261	Double_t neutralE = 0;
111318e8	262	Double_t maxCh = 0;
111318e8	263	Double_t maxNe = 0;
f5f3c8e9	264	Int_t gall = 0;
f5f3c8e9	265	Int_t gemc = 0;
111318e8	266	Int_t ncharged = 0;
111318e8	267	Int_t nneutral = 0;
02c7e943	268	Double_t mcpt = 0;
b65fac7a	269
f5f3c8e9	270	jet->SetNumberOfTracks(constituents.size());
f5f3c8e9	271	jet->SetNumberOfClusters(constituents.size());
b65fac7a	272	for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff	273	Int_t uid = constituents[ic].user_index();
b65fac7a	274
b65fac7a	275	if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
f5f3c8e9	276	++gall;
b65fac7a	277	Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
f5f3c8e9	278	Double_t geta = constituents[ic].eta();
b65fac7a	279	if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
b65fac7a	280	(geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
f5f3c8e9	281	++gemc;
f5f3c8e9	282	continue;
02c7e943	283	} else if ((uid > 0) && tracks) { // track constituent
111318e8	284	Int_t tid = uid - 100;
b65fac7a	285	AliVParticle t = static_cast<AliVParticle>(tracks->At(tid));
d03084cd	286	if (t) {
111318e8	287	if (t->Charge() == 0) {
b65fac7a	288	neutralE += t->P();
111318e8	289	++nneutral;
	290	if (t->Pt() > maxNe)
	291	maxNe = t->Pt();
	292	} else {
	293	++ncharged;
	294	if (t->Pt() > maxCh)
	295	maxCh = t->Pt();
	296	}
	297
	298	if (t->InheritsFrom("AliMCParticle") \|\| t->GetLabel() == 100) // MC particle
02c7e943	299	mcpt += t->Pt();
111318e8	300
b65fac7a	301	jet->AddTrackAt(tid, nt);
111318e8	302	++nt;
d03084cd	303	}
02c7e943	304	} else if (clus) { // cluster constituent
111318e8	305	Int_t cid = -uid - 100;
b65fac7a	306	AliVCluster c = static_cast<AliVCluster>(clus->At(cid));
	307	if (c) {
	308	TLorentzVector nP;
	309	c->GetMomentum(nP, vertex);
	310	neutralE += nP.P();
111318e8	311	if (nP.Pt() > maxNe)
	312	maxNe = nP.Pt();
	313
	314	if (c->Chi2() == 100) // MC particle
02c7e943	315	mcpt += nP.Pt();
111318e8	316
b65fac7a	317	jet->AddClusterAt(cid, nc);
111318e8	318	++nc;
111318e8	319	++nneutral;
b65fac7a	320	}
02c7e943	321	} else {
	322	AliError(Form("%s: No logical way to end up here.", GetName()));
	323	continue;
48d874ff	324	}
48d874ff	325	}
35789a2d	326	jet->SetNumberOfTracks(nt);
35789a2d	327	jet->SetNumberOfClusters(nc);
f5f3c8e9	328	jet->SortConstituents();
111318e8	329	jet->SetMaxNeutralPt(maxNe);
111318e8	330	jet->SetMaxChargedPt(maxCh);
b65fac7a	331	jet->SetNEF(neutralE / jet->E());
f5f3c8e9	332	jet->SetArea(fjw.GetJetArea(ij));
111318e8	333	jet->SetNumberOfCharged(ncharged);
111318e8	334	jet->SetNumberOfNeutrals(nneutral);
02c7e943	335	jet->SetMCPt(mcpt);
b65fac7a	336	if (gall > 0)
b65fac7a	337	jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9	338	else
f5f3c8e9	339	jet->SetAreaEmc(-1);
b65fac7a	340	if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65	341	(jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
	342	(jet->Eta() > geom->GetArm1EtaMin()) &&
	343	(jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9	344	jet->SetAxisInEmcal(kTRUE);
48d874ff	345	jetCount++;
	346	}
	347	}