[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetContainer.cxx

// $Id$
//
// Container with name, TClonesArray and cuts for jets
//
// Author: M. Verweij

#include <TClonesArray.h>

#include "AliEmcalJet.h"
#include "AliVEvent.h"
#include "AliLog.h"
#include "AliEMCALGeometry.h"
#include "AliParticleContainer.h"
#include "AliClusterContainer.h"
#include "AliLocalRhoParameter.h"

#include "AliJetContainer.h"

ClassImp(AliJetContainer)

//________________________________________________________________________
AliJetContainer::AliJetContainer():
  AliEmcalContainer("AliJetContainer"),
  fJetAcceptanceType(kUser),
  fJetRadius(0),
  fRhoName(),
  fLocalRhoName(),
  fFlavourSelection(0),
  fPtBiasJetTrack(0),
  fPtBiasJetClus(0),
  fJetPtCut(1),
  fJetAreaCut(-1),
  fAreaEmcCut(0),
  fJetMinEta(-0.9),
  fJetMaxEta(0.9),
  fJetMinPhi(-10),
  fJetMaxPhi(10),
  fMaxClusterPt(1000),
  fMaxTrackPt(100),
  fZLeadingEmcCut(10.),
  fZLeadingChCut(10.),
  fNEFMinCut(-10.),
  fNEFMaxCut(10.),
  fLeadingHadronType(0),
  fNLeadingJets(1),
  fJetBitMap(0),
  fJetTrigger(0),
  fParticleContainer(0),
  fClusterContainer(0),
  fRho(0),
  fLocalRho(0),
  fGeom(0),
  fRunNumber(0)
{
  // Default constructor.

  fClassName = "AliEmcalJet";
}

//________________________________________________________________________
AliJetContainer::AliJetContainer(const char *name):
  AliEmcalContainer(name),
  fJetAcceptanceType(kUser),
  fJetRadius(0),
  fRhoName(),
  fLocalRhoName(),
  fFlavourSelection(0),
  fPtBiasJetTrack(0),
  fPtBiasJetClus(0),
  fJetPtCut(1),
  fJetAreaCut(-1),
  fAreaEmcCut(0),
  fJetMinEta(-0.9),
  fJetMaxEta(0.9),
  fJetMinPhi(-10),
  fJetMaxPhi(10),
  fMaxClusterPt(1000),
  fMaxTrackPt(100),
  fZLeadingEmcCut(10.),
  fZLeadingChCut(10.),
  fNEFMinCut(-10.),
  fNEFMaxCut(10.),
  fLeadingHadronType(0),
  fNLeadingJets(1),
  fJetBitMap(0),
  fJetTrigger(0),
  fParticleContainer(0),
  fClusterContainer(0),
  fRho(0),
  fLocalRho(0),
  fGeom(0),
  fRunNumber(0)
{
  // Standard constructor.

  fClassName = "AliEmcalJet";
}

//________________________________________________________________________
void AliJetContainer::SetArray(AliVEvent *event) 
{
  // Set jet array

  AliEmcalContainer::SetArray(event);

  if(fJetAcceptanceType==kTPC) {
    AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
    SetJetEtaPhiTPC();
  }
  else if(fJetAcceptanceType==kEMCAL) {
    AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
    SetJetEtaPhiEMCAL();
 }
}


//________________________________________________________________________
void AliJetContainer::SetEMCALGeometry() {
  fGeom = AliEMCALGeometry::GetInstance();
  if (!fGeom) {
    AliError(Form("%s: Can not create geometry", GetName()));
    return;
  }
}

//________________________________________________________________________
void AliJetContainer::LoadRho(AliVEvent *event)
{
  // Load rho

  if (!fRhoName.IsNull() && !fRho) {
    fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
    if (!fRho) {
      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
      return;
    }
  }
}

//________________________________________________________________________
void AliJetContainer::LoadLocalRho(AliVEvent *event)
{
  // Load local rho

  if (!fLocalRhoName.IsNull() && !fLocalRho) {
    fLocalRho = dynamic_cast<AliLocalRhoParameter*>(event->FindListObject(fLocalRhoName));
    if (!fLocalRho) {
      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fLocalRhoName.Data()));
      return;
    }
  }
}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
{
  // Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho

  TString option(opt);
  option.ToLower();

  Int_t tempID = fCurrentID;

  AliEmcalJet *jetMax = GetNextAcceptJet(0);
  AliEmcalJet *jet = 0;

  if (option.Contains("rho")) {
    while ((jet = GetNextAcceptJet())) {
      if (jet->Pt()-jet->Area()*fRho->GetVal() > jetMax->Pt()-jetMax->Area()*fRho->GetVal()) jetMax = jet;
    }
  }
  else {
    while ((jet = GetNextAcceptJet())) {
      if (jet->Pt() > jetMax->Pt()) jetMax = jet;
    }
  }

  fCurrentID = tempID;

  return jetMax;
}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {

  //Get i^th jet in array

  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
  return jet;

}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {

  //Only return jet if is accepted

  AliEmcalJet *jet = GetJet(i);
  if(!AcceptJet(jet)) return 0;

  return jet;
}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {

  //Get particle with label lab in array
  
  Int_t i = GetIndexFromLabel(lab);
  return GetJet(i);
}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {

  //Get particle with label lab in array
  
  Int_t i = GetIndexFromLabel(lab);
  return GetAcceptJet(i);
}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {

  //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found

  if (i>=0) fCurrentID = i;

  const Int_t njets = GetNEntries();
  AliEmcalJet *jet = 0;
  while (fCurrentID < njets && !jet) { 
    jet = GetAcceptJet(fCurrentID);
    fCurrentID++;
  }

  return jet;
}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {

  //Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found

  if (i>=0) fCurrentID = i;

  const Int_t njets = GetNEntries();
  AliEmcalJet *jet = 0;
  while (fCurrentID < njets && !jet) { 
    jet = GetJet(fCurrentID);
    fCurrentID++;
  }

  return jet;
}

//________________________________________________________________________
Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
  AliEmcalJet *jet = GetJet(i);

  return jet->Pt() - fRho->GetVal()*jet->Area();
}

//________________________________________________________________________
Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
  AliEmcalJet *jet = GetJet(i);

  return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
}

//________________________________________________________________________
void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
{
  //Get momentum of the i^th jet in array

  AliEmcalJet *jet = GetJet(i);
  if(jet) jet->GetMom(mom);
}

//________________________________________________________________________
Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
{ 
  // Accept jet with a bias.

  if (fLeadingHadronType == 0) {
    if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
  }
  else if (fLeadingHadronType == 1) {
    if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
  }
  else {
    if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
  }

  return kTRUE;


}

//________________________________________________________________________
Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
{   

   // Return true if jet is accepted.

   if (!jet)
      return kFALSE;

   if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
      return kFALSE;

   if (jet->Pt() <= fJetPtCut) 
      return kFALSE;

   if (jet->Area() <= fJetAreaCut) 
      return kFALSE;

   if (jet->AreaEmc() < fAreaEmcCut)
      return kFALSE;
   
   if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
      return kFALSE;
   
   if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
      return kFALSE;

   if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut)
      return kFALSE;
   
   if (!AcceptBiasJet(jet))
      return kFALSE;
   
   if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
      return kFALSE;
   
   if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
      return kFALSE;
   
   Double_t jetPhi = jet->Phi();
   Double_t jetEta = jet->Eta();
   
   if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
      jetPhi -= TMath::Pi() * 2;
   
   return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
}

//________________________________________________________________________
Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
{
  if (fLeadingHadronType == 0)       // charged leading hadron
    return jet->MaxTrackPt();
  else if (fLeadingHadronType == 1)  // neutral leading hadron
    return jet->MaxClusterPt();
  else                               // charged or neutral
    return jet->MaxPartPt();
}

//________________________________________________________________________
void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
{
  Double_t maxClusterPt = 0;
  Double_t maxClusterEta = 0;
  Double_t maxClusterPhi = 0;
  
  Double_t maxTrackPt = 0;
  Double_t maxTrackEta = 0;
  Double_t maxTrackPhi = 0;
      
  if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
    AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
    if (cluster) {
      TLorentzVector nPart;
      cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
      
      maxClusterEta = nPart.Eta();
      maxClusterPhi = nPart.Phi();
      maxClusterPt = nPart.Pt();
    }
  }
      
  if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
    AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
    if (track) {
      maxTrackEta = track->Eta();
      maxTrackPhi = track->Phi();
      maxTrackPt = track->Pt();
    }
  }
      
  if (maxTrackPt > maxClusterPt) 
    mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
  else 
    mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
}

//________________________________________________________________________
Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
{

  if (fClusterContainer && fClusterContainer->GetArray()) {
    TLorentzVector mom;
    
    AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
    if (cluster) {
      cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
      
      return GetZ(jet,mom);
    }
    else
      return -1;
  }
  else
    return -1;
}

//________________________________________________________________________
Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
{

  if (fParticleContainer && fParticleContainer->GetArray() ) {
    TLorentzVector mom;
    
    AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
    if (track) {
      mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
      
      return GetZ(jet,mom);
    }
    else
      return -1;
  }
  else
    return -1;
}

//________________________________________________________________________
Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
{

  Double_t pJetSq = jet->Px()*jet->Px() + jet->Py()*jet->Py() + jet->Pz()*jet->Pz();

  if(pJetSq>1e-6)
    return (mom.Px()*jet->Px() + mom.Py()*jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
  else {
    AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
    return -1;
  }

}

//________________________________________________________________________
void AliJetContainer::SetJetEtaPhiEMCAL()
{
  //Set default cuts for full jets

  if(!fGeom) SetEMCALGeometry();
  if(fGeom) {
    SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);

    if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
      SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
    else
      SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);

  }
  else {
    AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
    SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
    SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
  }
}

//________________________________________________________________________
void AliJetContainer::SetJetEtaPhiTPC()
{
  //Set default cuts for charged jets

  SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
  SetJetPhiLimits(-10, 10);
}

//________________________________________________________________________
void AliJetContainer::ResetCuts() 
{
  // Reset cuts to default values

  fPtBiasJetTrack = 0;
  fPtBiasJetClus = 0;
  fJetPtCut = 1;
  fJetAreaCut = -1;
  fAreaEmcCut = 0;
  fJetMinEta = -0.9;
  fJetMaxEta = 0.9;
  fJetMinPhi = -10;
  fJetMaxPhi = 10;
  fMaxClusterPt = 1000;
  fMaxTrackPt = 100;
  fLeadingHadronType = 0;
  fZLeadingEmcCut = 10.;
  fZLeadingChCut = 10.;
}

//________________________________________________________________________
void AliJetContainer::SetClassName(const char *clname)
{
  // Set the class name

  TClass cls(clname);
  if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
  else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
}
Commit	Line	Data
9239b066	1	// $Id$
421b12a3	2	//
9239b066	3	// Container with name, TClonesArray and cuts for jets
421b12a3	4	//
	5	// Author: M. Verweij
	6
421b12a3	7	#include <TClonesArray.h>
6421eeb0	8
421b12a3	9	#include "AliEmcalJet.h"
	10	#include "AliVEvent.h"
	11	#include "AliLog.h"
421b12a3	12	#include "AliEMCALGeometry.h"
6421eeb0	13	#include "AliParticleContainer.h"
6421eeb0	14	#include "AliClusterContainer.h"
9e5eee5d	15	#include "AliLocalRhoParameter.h"
421b12a3	16
	17	#include "AliJetContainer.h"
	18
	19	ClassImp(AliJetContainer)
	20
	21	//________________________________________________________________________
	22	AliJetContainer::AliJetContainer():
	23	AliEmcalContainer("AliJetContainer"),
	24	fJetAcceptanceType(kUser),
	25	fJetRadius(0),
	26	fRhoName(),
9e5eee5d	27	fLocalRhoName(),
3c9775d9	28	fFlavourSelection(0),
421b12a3	29	fPtBiasJetTrack(0),
	30	fPtBiasJetClus(0),
	31	fJetPtCut(1),
	32	fJetAreaCut(-1),
	33	fAreaEmcCut(0),
	34	fJetMinEta(-0.9),
	35	fJetMaxEta(0.9),
	36	fJetMinPhi(-10),
	37	fJetMaxPhi(10),
	38	fMaxClusterPt(1000),
	39	fMaxTrackPt(100),
413ab3a3	40	fZLeadingEmcCut(10.),
413ab3a3	41	fZLeadingChCut(10.),
9e5eee5d	42	fNEFMinCut(-10.),
9e5eee5d	43	fNEFMaxCut(10.),
421b12a3	44	fLeadingHadronType(0),
	45	fNLeadingJets(1),
	46	fJetBitMap(0),
6421eeb0	47	fJetTrigger(0),
	48	fParticleContainer(0),
	49	fClusterContainer(0),
421b12a3	50	fRho(0),
9e5eee5d	51	fLocalRho(0),
56f370b0	52	fGeom(0),
56f370b0	53	fRunNumber(0)
421b12a3	54	{
	55	// Default constructor.
	56
b6f970ad	57	fClassName = "AliEmcalJet";
421b12a3	58	}
	59
	60	//________________________________________________________________________
	61	AliJetContainer::AliJetContainer(const char *name):
	62	AliEmcalContainer(name),
	63	fJetAcceptanceType(kUser),
	64	fJetRadius(0),
	65	fRhoName(),
9e5eee5d	66	fLocalRhoName(),
3c9775d9	67	fFlavourSelection(0),
421b12a3	68	fPtBiasJetTrack(0),
	69	fPtBiasJetClus(0),
	70	fJetPtCut(1),
	71	fJetAreaCut(-1),
	72	fAreaEmcCut(0),
	73	fJetMinEta(-0.9),
	74	fJetMaxEta(0.9),
	75	fJetMinPhi(-10),
	76	fJetMaxPhi(10),
	77	fMaxClusterPt(1000),
	78	fMaxTrackPt(100),
413ab3a3	79	fZLeadingEmcCut(10.),
413ab3a3	80	fZLeadingChCut(10.),
9e5eee5d	81	fNEFMinCut(-10.),
9e5eee5d	82	fNEFMaxCut(10.),
421b12a3	83	fLeadingHadronType(0),
	84	fNLeadingJets(1),
	85	fJetBitMap(0),
6421eeb0	86	fJetTrigger(0),
	87	fParticleContainer(0),
	88	fClusterContainer(0),
421b12a3	89	fRho(0),
9e5eee5d	90	fLocalRho(0),
56f370b0	91	fGeom(0),
56f370b0	92	fRunNumber(0)
421b12a3	93	{
	94	// Standard constructor.
	95
b6f970ad	96	fClassName = "AliEmcalJet";
421b12a3	97	}
421b12a3	98
421b12a3	99	//________________________________________________________________________
b6f970ad	100	void AliJetContainer::SetArray(AliVEvent *event)
421b12a3	101	{
	102	// Set jet array
	103
b6f970ad	104	AliEmcalContainer::SetArray(event);
421b12a3	105
56f370b0	106	if(fJetAcceptanceType==kTPC) {
56f370b0	107	AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
421b12a3	108	SetJetEtaPhiTPC();
56f370b0	109	}
	110	else if(fJetAcceptanceType==kEMCAL) {
	111	AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
421b12a3	112	SetJetEtaPhiEMCAL();
56f370b0	113	}
421b12a3	114	}
421b12a3	115
56f370b0	116
421b12a3	117	//________________________________________________________________________
	118	void AliJetContainer::SetEMCALGeometry() {
	119	fGeom = AliEMCALGeometry::GetInstance();
	120	if (!fGeom) {
	121	AliError(Form("%s: Can not create geometry", GetName()));
	122	return;
	123	}
	124	}
	125
	126	//________________________________________________________________________
	127	void AliJetContainer::LoadRho(AliVEvent *event)
	128	{
6421eeb0	129	// Load rho
421b12a3	130
	131	if (!fRhoName.IsNull() && !fRho) {
	132	fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
	133	if (!fRho) {
	134	AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
	135	return;
	136	}
	137	}
	138	}
	139
9e5eee5d	140	//________________________________________________________________________
	141	void AliJetContainer::LoadLocalRho(AliVEvent *event)
	142	{
	143	// Load local rho
	144
	145	if (!fLocalRhoName.IsNull() && !fLocalRho) {
	146	fLocalRho = dynamic_cast<AliLocalRhoParameter*>(event->FindListObject(fLocalRhoName));
	147	if (!fLocalRho) {
	148	AliError(Form("%s: Could not retrieve rho %s!", GetName(), fLocalRhoName.Data()));
	149	return;
	150	}
	151	}
	152	}
	153
6421eeb0	154	//________________________________________________________________________
b6f970ad	155	AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
6421eeb0	156	{
	157	// Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
	158
	159	TString option(opt);
	160	option.ToLower();
	161
b6f970ad	162	Int_t tempID = fCurrentID;
b6f970ad	163
6421eeb0	164	AliEmcalJet *jetMax = GetNextAcceptJet(0);
	165	AliEmcalJet *jet = 0;
	166
	167	if (option.Contains("rho")) {
	168	while ((jet = GetNextAcceptJet())) {
	169	if (jet->Pt()-jet->Area()fRho->GetVal() > jetMax->Pt()-jetMax->Area()fRho->GetVal()) jetMax = jet;
	170	}
	171	}
	172	else {
	173	while ((jet = GetNextAcceptJet())) {
	174	if (jet->Pt() > jetMax->Pt()) jetMax = jet;
	175	}
	176	}
	177
b6f970ad	178	fCurrentID = tempID;
b6f970ad	179
6421eeb0	180	return jetMax;
	181	}
	182
421b12a3	183	//________________________________________________________________________
	184	AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
	185
	186	//Get i^th jet in array
	187
	188	if(i<0 \|\| i>fClArray->GetEntriesFast()) return 0;
	189	AliEmcalJet jet = static_cast<AliEmcalJet>(fClArray->At(i));
	190	return jet;
	191
	192	}
	193
421b12a3	194	//________________________________________________________________________
	195	AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
	196
	197	//Only return jet if is accepted
	198
	199	AliEmcalJet *jet = GetJet(i);
	200	if(!AcceptJet(jet)) return 0;
	201
	202	return jet;
	203	}
	204
7cd832c7	205	//________________________________________________________________________
	206	AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
	207
	208	//Get particle with label lab in array
	209
	210	Int_t i = GetIndexFromLabel(lab);
	211	return GetJet(i);
	212	}
	213
	214	//________________________________________________________________________
	215	AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {
	216
	217	//Get particle with label lab in array
	218
	219	Int_t i = GetIndexFromLabel(lab);
	220	return GetAcceptJet(i);
	221	}
	222
6421eeb0	223	//________________________________________________________________________
b6f970ad	224	AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
6421eeb0	225
	226	//Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
	227
b6f970ad	228	if (i>=0) fCurrentID = i;
6421eeb0	229
	230	const Int_t njets = GetNEntries();
	231	AliEmcalJet *jet = 0;
b6f970ad	232	while (fCurrentID < njets && !jet) {
	233	jet = GetAcceptJet(fCurrentID);
	234	fCurrentID++;
6421eeb0	235	}
	236
	237	return jet;
	238	}
	239
	240	//________________________________________________________________________
7cd832c7	241	AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
	242
	243	//Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
	244
	245	if (i>=0) fCurrentID = i;
	246
	247	const Int_t njets = GetNEntries();
	248	AliEmcalJet *jet = 0;
	249	while (fCurrentID < njets && !jet) {
	250	jet = GetJet(fCurrentID);
	251	fCurrentID++;
	252	}
	253
	254	return jet;
	255	}
	256
	257	//________________________________________________________________________
	258	Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
	259	AliEmcalJet *jet = GetJet(i);
	260
	261	return jet->Pt() - fRho->GetVal()*jet->Area();
	262	}
	263
9e5eee5d	264	//________________________________________________________________________
	265	Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
	266	AliEmcalJet *jet = GetJet(i);
	267
	268	return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
	269	}
	270
7cd832c7	271	//________________________________________________________________________
6421eeb0	272	void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
	273	{
	274	//Get momentum of the i^th jet in array
	275
	276	AliEmcalJet *jet = GetJet(i);
	277	if(jet) jet->GetMom(mom);
	278	}
	279
421b12a3	280	//________________________________________________________________________
	281	Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
	282	{
	283	// Accept jet with a bias.
	284
	285	if (fLeadingHadronType == 0) {
	286	if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
	287	}
	288	else if (fLeadingHadronType == 1) {
	289	if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
	290	}
	291	else {
	292	if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
	293	}
	294
	295	return kTRUE;
	296
	297
	298	}
	299
	300	//________________________________________________________________________
	301	Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
	302	{
421b12a3	303
3c9775d9	304	// Return true if jet is accepted.
	305
	306	if (!jet)
	307	return kFALSE;
	308
	309	if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
	310	return kFALSE;
	311
	312	if (jet->Pt() <= fJetPtCut)
	313	return kFALSE;
	314
	315	if (jet->Area() <= fJetAreaCut)
	316	return kFALSE;
	317
	318	if (jet->AreaEmc() < fAreaEmcCut)
	319	return kFALSE;
	320
	321	if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
	322	return kFALSE;
	323
	324	if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
	325	return kFALSE;
	326
	327	if (jet->NEF() < fNEFMinCut \|\| jet->NEF() > fNEFMaxCut)
	328	return kFALSE;
	329
	330	if (!AcceptBiasJet(jet))
	331	return kFALSE;
	332
	333	if (jet->MaxTrackPt() > fMaxTrackPt \|\| jet->MaxClusterPt() > fMaxClusterPt)
	334	return kFALSE;
	335
	336	if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
	337	return kFALSE;
	338
	339	Double_t jetPhi = jet->Phi();
	340	Double_t jetEta = jet->Eta();
	341
	342	if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
	343	jetPhi -= TMath::Pi() * 2;
	344
	345	return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
421b12a3	346	}
	347
	348	//________________________________________________________________________
	349	Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
	350	{
	351	if (fLeadingHadronType == 0) // charged leading hadron
	352	return jet->MaxTrackPt();
	353	else if (fLeadingHadronType == 1) // neutral leading hadron
	354	return jet->MaxClusterPt();
	355	else // charged or neutral
	356	return jet->MaxPartPt();
	357	}
	358
6421eeb0	359	//________________________________________________________________________
	360	void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
	361	{
	362	Double_t maxClusterPt = 0;
	363	Double_t maxClusterEta = 0;
	364	Double_t maxClusterPhi = 0;
	365
	366	Double_t maxTrackPt = 0;
	367	Double_t maxTrackEta = 0;
	368	Double_t maxTrackPhi = 0;
	369
	370	if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 \|\| fLeadingHadronType == 2)) {
	371	AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
	372	if (cluster) {
	373	TLorentzVector nPart;
	374	cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
	375
	376	maxClusterEta = nPart.Eta();
	377	maxClusterPhi = nPart.Phi();
	378	maxClusterPt = nPart.Pt();
	379	}
	380	}
	381
	382	if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 \|\| fLeadingHadronType == 2)) {
	383	AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
	384	if (track) {
	385	maxTrackEta = track->Eta();
	386	maxTrackPhi = track->Phi();
	387	maxTrackPt = track->Pt();
	388	}
	389	}
	390
	391	if (maxTrackPt > maxClusterPt)
	392	mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
	393	else
	394	mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
	395	}
421b12a3	396
413ab3a3	397	//________________________________________________________________________
	398	Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
	399	{
	400
ad5d2c11	401	if (fClusterContainer && fClusterContainer->GetArray()) {
413ab3a3	402	TLorentzVector mom;
	403
	404	AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
	405	if (cluster) {
	406	cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
	407
	408	return GetZ(jet,mom);
	409	}
	410	else
	411	return -1;
	412	}
	413	else
	414	return -1;
	415	}
	416
	417	//________________________________________________________________________
	418	Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
	419	{
	420
	421	if (fParticleContainer && fParticleContainer->GetArray() ) {
	422	TLorentzVector mom;
	423
	424	AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
	425	if (track) {
	426	mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
	427
	428	return GetZ(jet,mom);
	429	}
	430	else
	431	return -1;
	432	}
	433	else
	434	return -1;
	435	}
	436
	437	//________________________________________________________________________
	438	Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
	439	{
	440
	441	Double_t pJetSq = jet->Px()jet->Px() + jet->Py()jet->Py() + jet->Pz()*jet->Pz();
	442
ad5d2c11	443	if(pJetSq>1e-6)
413ab3a3	444	return (mom.Px()jet->Px() + mom.Py()jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
	445	else {
	446	AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
	447	return -1;
	448	}
	449
	450	}
	451
421b12a3	452	//________________________________________________________________________
	453	void AliJetContainer::SetJetEtaPhiEMCAL()
	454	{
	455	//Set default cuts for full jets
	456
	457	if(!fGeom) SetEMCALGeometry();
	458	if(fGeom) {
	459	SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
56f370b0	460
	461	if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
	462	SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
	463	else
	464	SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
	465
421b12a3	466	}
421b12a3	467	else {
56f370b0	468	AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
421b12a3	469	SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
254370e1	470	SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
421b12a3	471	}
421b12a3	472	}
	473
	474	//________________________________________________________________________
	475	void AliJetContainer::SetJetEtaPhiTPC()
	476	{
56f370b0	477	//Set default cuts for charged jets
421b12a3	478
56f370b0	479	SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
421b12a3	480	SetJetPhiLimits(-10, 10);
421b12a3	481	}
	482
	483	//________________________________________________________________________
	484	void AliJetContainer::ResetCuts()
	485	{
	486	// Reset cuts to default values
	487
	488	fPtBiasJetTrack = 0;
	489	fPtBiasJetClus = 0;
	490	fJetPtCut = 1;
	491	fJetAreaCut = -1;
	492	fAreaEmcCut = 0;
	493	fJetMinEta = -0.9;
	494	fJetMaxEta = 0.9;
	495	fJetMinPhi = -10;
	496	fJetMaxPhi = 10;
	497	fMaxClusterPt = 1000;
	498	fMaxTrackPt = 100;
	499	fLeadingHadronType = 0;
413ab3a3	500	fZLeadingEmcCut = 10.;
413ab3a3	501	fZLeadingChCut = 10.;
421b12a3	502	}
b6f970ad	503
	504	//________________________________________________________________________
	505	void AliJetContainer::SetClassName(const char *clname)
	506	{
	507	// Set the class name
	508
	509	TClass cls(clname);
	510	if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
	511	else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
	512	}