[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 "AliPythiaInfo.h"

#include "AliJetContainer.h"

ClassImp(AliJetContainer)

//________________________________________________________________________
AliJetContainer::AliJetContainer():
  AliEmcalContainer("AliJetContainer"),
  fJetAcceptanceType(kUser),
  fJetRadius(0),
  fRhoName(),
  fLocalRhoName(),
  fRhoMassName(),
  fPythiaInfoName(),
  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),
  fTagStatus(-1),
  fParticleContainer(0),
  fClusterContainer(0),
  fRho(0),
  fLocalRho(0),
  fRhoMass(0),
  fPythiaInfo(0),
  fGeom(0),
  fRunNumber(0)
{
  // Default constructor.

  fClassName = "AliEmcalJet";
}

//________________________________________________________________________
AliJetContainer::AliJetContainer(const char *name):
  AliEmcalContainer(name),
  fJetAcceptanceType(kUser),
  fJetRadius(0),
  fRhoName(),
  fLocalRhoName(),
  fRhoMassName(),
  fPythiaInfoName(),
  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),
  fTagStatus(-1),
  fParticleContainer(0),
  fClusterContainer(0),
  fRho(0),
  fLocalRho(0),
  fRhoMass(0),
  fPythiaInfo(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;
    }
  }
}

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

  if (!fRhoMassName.IsNull() && !fRhoMass) {
    fRhoMass = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoMassName));
    if (!fRhoMass) {
      AliError(Form("%s: Could not retrieve rho_mass %s!", GetName(), fRhoMassName.Data()));
      return;
    }
  }
}
//________________________________________________________________________
void AliJetContainer::LoadPythiaInfo(AliVEvent *event)
{
    // Load parton info
    
    if (!fPythiaInfoName.IsNull() && !fPythiaInfo) {
        fPythiaInfo = dynamic_cast<AliPythiaInfo*>(event->FindListObject(fPythiaInfoName));
        if (!fPythiaInfo) {
           AliError(Form("%s: Could not retrieve parton infos! %s!", GetName(), fPythiaInfoName.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()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) )
	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) {

  //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) {

  //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(const AliEmcalJet *jet)
{ 
  // 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(const AliEmcalJet *jet)
{   
   // Return true if jet is accepted.

  fRejectionReason = 0;

  if (!jet) {
    AliDebug(11,"No jet found");
    fRejectionReason |= kNullObject;
    return kFALSE;
  }

  if (jet->Pt() <= fJetPtCut) {
    AliDebug(11,Form("Cut rejecting jet: JetPtCut %.1f",fJetPtCut));
    fRejectionReason |= kPtCut;
    return kFALSE;
  }

  Double_t jetPhi = jet->Phi();
  Double_t jetEta = jet->Eta();
   
  // if limits are given in (-pi, pi) range
  if (fJetMinPhi < 0) jetPhi -= TMath::Pi() * 2;
   
  if (jetEta < fJetMinEta || jetEta > fJetMaxEta || jetPhi < fJetMinPhi || jetPhi > fJetMaxPhi) {
    AliDebug(11,"Cut rejecting jet: Acceptance");
    fRejectionReason |= kAcceptanceCut;
    return kFALSE;
  }

  if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap) {
    AliDebug(11,"Cut rejecting jet: Bit map");
    fRejectionReason |= kBitMapCut;
    return kFALSE;
  }

  if (jet->Area() <= fJetAreaCut)  {
    AliDebug(11,"Cut rejecting jet: Area");
    fRejectionReason |= kAreaCut;
    return kFALSE;
  }

  if (jet->AreaEmc() < fAreaEmcCut) {
    AliDebug(11,"Cut rejecting jet: AreaEmc");
    fRejectionReason |= kAreaEmcCut;
    return kFALSE;
  }
   
  if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut) {
    AliDebug(11,"Cut rejecting jet: ZLeading");
    fRejectionReason |= kZLeadingChCut;
    return kFALSE;
  }
   
  if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut) {
    AliDebug(11,"Cut rejecting jet: ZLeadEmc");
    fRejectionReason |= kZLeadingEmcCut;
    return kFALSE;
  }

  if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut) {
    AliDebug(11,"Cut rejecting jet: NEF");
    fRejectionReason |= kNEFCut;
    return kFALSE;
  }
   
  if (!AcceptBiasJet(jet)) {
    AliDebug(11,"Cut rejecting jet: Bias");
    fRejectionReason |= kMinLeadPtCut;
    return kFALSE;
  }

  if (jet->MaxTrackPt() > fMaxTrackPt) {
    AliDebug(11,"Cut rejecting jet: MaxTrackPt");
    fRejectionReason |= kMaxTrackPtCut;
    return kFALSE;

  }

  if (jet->MaxClusterPt() > fMaxClusterPt) {
    AliDebug(11,"Cut rejecting jet: MaxClusPt");
    fRejectionReason |= kMaxClusterPtCut;
    return kFALSE;
  }
   
  if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection)) {
    AliDebug(11,"Cut rejecting jet: Flavour");
    fRejectionReason |= kFlavourCut;
    return kFALSE;
  }
   
  if(fTagStatus>-1 && jet->GetTagStatus()!=fTagStatus) {
    AliDebug(11,"Cut rejecting jet: tag status");
    fRejectionReason |= kTagStatus;
    return kFALSE;
  }

  return kTRUE;
}

//________________________________________________________________________
Double_t AliJetContainer::GetLeadingHadronPt(const 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, const 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(const 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(const 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(const 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.405+fJetRadius,3.135-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::PrintCuts() 
{
  // Reset cuts to default values

  Printf("PtBiasJetTrack: %f",fPtBiasJetTrack);
  Printf("PtBiasJetClus: %f",fPtBiasJetClus);
  Printf("JetPtCut: %f", fJetPtCut);
  Printf("JetAreaCut: %f",fJetAreaCut);
  Printf("AreaEmcCut: %f",fAreaEmcCut);
  Printf("JetMinEta: %f", fJetMinEta);
  Printf("JetMaxEta: %f", fJetMaxEta);
  Printf("JetMinPhi: %f", fJetMinPhi);
  Printf("JetMaxPhi: %f", fJetMaxPhi);
  Printf("MaxClusterPt: %f",fMaxClusterPt);
  Printf("MaxTrackPt: %f",fMaxTrackPt);
  Printf("LeadingHadronType: %d",fLeadingHadronType);
  Printf("ZLeadingEmcCut: %f",fZLeadingEmcCut);
  Printf("ZLeadingChCut: %f",fZLeadingChCut);

}

//________________________________________________________________________
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));
}

//________________________________________________________________________
Double_t AliJetContainer::GetFractionSharedPt(const AliEmcalJet *jet1) const
{
  //
  // Get fraction of shared pT between matched jets
  // Uses ClosestJet() jet pT as baseline: fraction = \Sum_{const,jet1} pT,const,i / pT,jet,closest
  // Only works if tracks array of both jets is the same
  //

  AliEmcalJet *jet2 = jet1->ClosestJet();
  if(!jet2) return -1;

  Double_t fraction = 0.;
  Double_t jetPt2 = jet2->Pt();
 
  if(jetPt2>0) {
    Double_t sumPt = 0.;
    AliVParticle *vpf = 0x0;
    Int_t iFound = 0;
    for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
      Int_t idx = (Int_t)jet2->TrackAt(icc);
      iFound = 0;
      for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
	if(idx == jet1->TrackAt(icf) && iFound==0 ) {
	  iFound=1;
	  vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fParticleContainer->GetArray()));
	  if(vpf) sumPt += vpf->Pt();
	  continue;
	}
      }
    }
    fraction = sumPt/jetPt2;
  } else 
    fraction = -1;
  
  return fraction;
}
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"
b929303d	16	#include "AliPythiaInfo.h"
421b12a3	17
	18	#include "AliJetContainer.h"
	19
	20	ClassImp(AliJetContainer)
	21
	22	//________________________________________________________________________
	23	AliJetContainer::AliJetContainer():
	24	AliEmcalContainer("AliJetContainer"),
	25	fJetAcceptanceType(kUser),
	26	fJetRadius(0),
	27	fRhoName(),
9e5eee5d	28	fLocalRhoName(),
4d3b366f	29	fRhoMassName(),
b929303d	30	fPythiaInfoName(),
3c9775d9	31	fFlavourSelection(0),
421b12a3	32	fPtBiasJetTrack(0),
	33	fPtBiasJetClus(0),
	34	fJetPtCut(1),
	35	fJetAreaCut(-1),
	36	fAreaEmcCut(0),
	37	fJetMinEta(-0.9),
	38	fJetMaxEta(0.9),
	39	fJetMinPhi(-10),
	40	fJetMaxPhi(10),
	41	fMaxClusterPt(1000),
	42	fMaxTrackPt(100),
413ab3a3	43	fZLeadingEmcCut(10.),
413ab3a3	44	fZLeadingChCut(10.),
9e5eee5d	45	fNEFMinCut(-10.),
9e5eee5d	46	fNEFMaxCut(10.),
421b12a3	47	fLeadingHadronType(0),
	48	fNLeadingJets(1),
	49	fJetBitMap(0),
6421eeb0	50	fJetTrigger(0),
7a615dc4	51	fTagStatus(-1),
6421eeb0	52	fParticleContainer(0),
6421eeb0	53	fClusterContainer(0),
421b12a3	54	fRho(0),
9e5eee5d	55	fLocalRho(0),
4d3b366f	56	fRhoMass(0),
b929303d	57	fPythiaInfo(0),
56f370b0	58	fGeom(0),
56f370b0	59	fRunNumber(0)
421b12a3	60	{
	61	// Default constructor.
	62
b6f970ad	63	fClassName = "AliEmcalJet";
421b12a3	64	}
	65
	66	//________________________________________________________________________
	67	AliJetContainer::AliJetContainer(const char *name):
	68	AliEmcalContainer(name),
	69	fJetAcceptanceType(kUser),
	70	fJetRadius(0),
	71	fRhoName(),
9e5eee5d	72	fLocalRhoName(),
4d3b366f	73	fRhoMassName(),
b929303d	74	fPythiaInfoName(),
3c9775d9	75	fFlavourSelection(0),
421b12a3	76	fPtBiasJetTrack(0),
	77	fPtBiasJetClus(0),
	78	fJetPtCut(1),
	79	fJetAreaCut(-1),
	80	fAreaEmcCut(0),
	81	fJetMinEta(-0.9),
	82	fJetMaxEta(0.9),
	83	fJetMinPhi(-10),
	84	fJetMaxPhi(10),
	85	fMaxClusterPt(1000),
	86	fMaxTrackPt(100),
413ab3a3	87	fZLeadingEmcCut(10.),
413ab3a3	88	fZLeadingChCut(10.),
9e5eee5d	89	fNEFMinCut(-10.),
9e5eee5d	90	fNEFMaxCut(10.),
421b12a3	91	fLeadingHadronType(0),
	92	fNLeadingJets(1),
	93	fJetBitMap(0),
6421eeb0	94	fJetTrigger(0),
7a615dc4	95	fTagStatus(-1),
6421eeb0	96	fParticleContainer(0),
6421eeb0	97	fClusterContainer(0),
421b12a3	98	fRho(0),
9e5eee5d	99	fLocalRho(0),
4d3b366f	100	fRhoMass(0),
b929303d	101	fPythiaInfo(0),
56f370b0	102	fGeom(0),
56f370b0	103	fRunNumber(0)
421b12a3	104	{
	105	// Standard constructor.
	106
b6f970ad	107	fClassName = "AliEmcalJet";
421b12a3	108	}
421b12a3	109
421b12a3	110	//________________________________________________________________________
b6f970ad	111	void AliJetContainer::SetArray(AliVEvent *event)
421b12a3	112	{
	113	// Set jet array
	114
b6f970ad	115	AliEmcalContainer::SetArray(event);
421b12a3	116
56f370b0	117	if(fJetAcceptanceType==kTPC) {
56f370b0	118	AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
421b12a3	119	SetJetEtaPhiTPC();
56f370b0	120	}
	121	else if(fJetAcceptanceType==kEMCAL) {
	122	AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
421b12a3	123	SetJetEtaPhiEMCAL();
56f370b0	124	}
421b12a3	125	}
421b12a3	126
56f370b0	127
421b12a3	128	//________________________________________________________________________
	129	void AliJetContainer::SetEMCALGeometry() {
	130	fGeom = AliEMCALGeometry::GetInstance();
	131	if (!fGeom) {
	132	AliError(Form("%s: Can not create geometry", GetName()));
	133	return;
	134	}
	135	}
	136
	137	//________________________________________________________________________
	138	void AliJetContainer::LoadRho(AliVEvent *event)
	139	{
6421eeb0	140	// Load rho
421b12a3	141
	142	if (!fRhoName.IsNull() && !fRho) {
	143	fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
	144	if (!fRho) {
	145	AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
	146	return;
	147	}
	148	}
	149	}
	150
9e5eee5d	151	//________________________________________________________________________
	152	void AliJetContainer::LoadLocalRho(AliVEvent *event)
	153	{
	154	// Load local rho
	155
	156	if (!fLocalRhoName.IsNull() && !fLocalRho) {
	157	fLocalRho = dynamic_cast<AliLocalRhoParameter*>(event->FindListObject(fLocalRhoName));
	158	if (!fLocalRho) {
	159	AliError(Form("%s: Could not retrieve rho %s!", GetName(), fLocalRhoName.Data()));
	160	return;
	161	}
	162	}
	163	}
	164
4d3b366f	165	//________________________________________________________________________
	166	void AliJetContainer::LoadRhoMass(AliVEvent *event)
	167	{
	168	// Load rho
	169
	170	if (!fRhoMassName.IsNull() && !fRhoMass) {
	171	fRhoMass = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoMassName));
	172	if (!fRhoMass) {
	173	AliError(Form("%s: Could not retrieve rho_mass %s!", GetName(), fRhoMassName.Data()));
	174	return;
	175	}
	176	}
	177	}
cfa7b4df	178	//________________________________________________________________________
b929303d	179	void AliJetContainer::LoadPythiaInfo(AliVEvent *event)
cfa7b4df	180	{
	181	// Load parton info
	182
b929303d	183	if (!fPythiaInfoName.IsNull() && !fPythiaInfo) {
	184	fPythiaInfo = dynamic_cast<AliPythiaInfo*>(event->FindListObject(fPythiaInfoName));
	185	if (!fPythiaInfo) {
	186	AliError(Form("%s: Could not retrieve parton infos! %s!", GetName(), fPythiaInfoName.Data()));
	187	return;
cfa7b4df	188	}
	189	}
	190	}
	191
4d3b366f	192
4d3b366f	193
6421eeb0	194	//________________________________________________________________________
b6f970ad	195	AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
6421eeb0	196	{
	197	// Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
	198
	199	TString option(opt);
	200	option.ToLower();
	201
b6f970ad	202	Int_t tempID = fCurrentID;
b6f970ad	203
6421eeb0	204	AliEmcalJet *jetMax = GetNextAcceptJet(0);
	205	AliEmcalJet *jet = 0;
	206
	207	if (option.Contains("rho")) {
	208	while ((jet = GetNextAcceptJet())) {
9f4cfa2c MV	209	if ( (jet->Pt()-jet->Area()GetRhoVal()) > (jetMax->Pt()-jetMax->Area()GetRhoVal()) )
9f4cfa2c MV	210	jetMax = jet;
6421eeb0	211	}
	212	}
	213	else {
	214	while ((jet = GetNextAcceptJet())) {
	215	if (jet->Pt() > jetMax->Pt()) jetMax = jet;
	216	}
	217	}
	218
b6f970ad	219	fCurrentID = tempID;
b6f970ad	220
6421eeb0	221	return jetMax;
	222	}
	223
421b12a3	224	//________________________________________________________________________
	225	AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
	226
	227	//Get i^th jet in array
	228
	229	if(i<0 \|\| i>fClArray->GetEntriesFast()) return 0;
	230	AliEmcalJet jet = static_cast<AliEmcalJet>(fClArray->At(i));
	231	return jet;
	232
	233	}
	234
421b12a3	235	//________________________________________________________________________
6504c17f	236	AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) {
421b12a3	237
	238	//Only return jet if is accepted
	239
	240	AliEmcalJet *jet = GetJet(i);
	241	if(!AcceptJet(jet)) return 0;
	242
	243	return jet;
	244	}
	245
7cd832c7	246	//________________________________________________________________________
	247	AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
	248
	249	//Get particle with label lab in array
	250
	251	Int_t i = GetIndexFromLabel(lab);
	252	return GetJet(i);
	253	}
	254
	255	//________________________________________________________________________
6504c17f	256	AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) {
7cd832c7	257
	258	//Get particle with label lab in array
	259
	260	Int_t i = GetIndexFromLabel(lab);
	261	return GetAcceptJet(i);
	262	}
	263
6421eeb0	264	//________________________________________________________________________
b6f970ad	265	AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
6421eeb0	266
	267	//Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
	268
b6f970ad	269	if (i>=0) fCurrentID = i;
6421eeb0	270
	271	const Int_t njets = GetNEntries();
	272	AliEmcalJet *jet = 0;
b6f970ad	273	while (fCurrentID < njets && !jet) {
	274	jet = GetAcceptJet(fCurrentID);
	275	fCurrentID++;
6421eeb0	276	}
	277
	278	return jet;
	279	}
	280
	281	//________________________________________________________________________
7cd832c7	282	AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
	283
	284	//Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
	285
	286	if (i>=0) fCurrentID = i;
	287
	288	const Int_t njets = GetNEntries();
	289	AliEmcalJet *jet = 0;
	290	while (fCurrentID < njets && !jet) {
	291	jet = GetJet(fCurrentID);
	292	fCurrentID++;
	293	}
	294
	295	return jet;
	296	}
	297
	298	//________________________________________________________________________
	299	Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
	300	AliEmcalJet *jet = GetJet(i);
	301
	302	return jet->Pt() - fRho->GetVal()*jet->Area();
	303	}
	304
9e5eee5d	305	//________________________________________________________________________
	306	Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
	307	AliEmcalJet *jet = GetJet(i);
	308
	309	return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
	310	}
	311
7cd832c7	312	//________________________________________________________________________
6421eeb0	313	void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
	314	{
	315	//Get momentum of the i^th jet in array
	316
	317	AliEmcalJet *jet = GetJet(i);
	318	if(jet) jet->GetMom(mom);
	319	}
	320
421b12a3	321	//________________________________________________________________________
983c25fc	322	Bool_t AliJetContainer::AcceptBiasJet(const AliEmcalJet *jet)
421b12a3	323	{
	324	// Accept jet with a bias.
	325
	326	if (fLeadingHadronType == 0) {
	327	if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
	328	}
	329	else if (fLeadingHadronType == 1) {
	330	if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
	331	}
	332	else {
	333	if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
	334	}
	335
	336	return kTRUE;
421b12a3	337	}
	338
	339	//________________________________________________________________________
983c25fc	340	Bool_t AliJetContainer::AcceptJet(const AliEmcalJet *jet)
421b12a3	341	{
3c9775d9	342	// Return true if jet is accepted.
3c9775d9	343
6504c17f	344	fRejectionReason = 0;
6504c17f	345
6c3dc2d6	346	if (!jet) {
6c3dc2d6	347	AliDebug(11,"No jet found");
6504c17f	348	fRejectionReason \|= kNullObject;
6c3dc2d6	349	return kFALSE;
6c3dc2d6	350	}
3c9775d9	351
6c3dc2d6	352	if (jet->Pt() <= fJetPtCut) {
6504c17f	353	AliDebug(11,Form("Cut rejecting jet: JetPtCut %.1f",fJetPtCut));
	354	fRejectionReason \|= kPtCut;
	355	return kFALSE;
	356	}
	357
	358	Double_t jetPhi = jet->Phi();
	359	Double_t jetEta = jet->Eta();
	360
	361	// if limits are given in (-pi, pi) range
	362	if (fJetMinPhi < 0) jetPhi -= TMath::Pi() * 2;
	363
	364	if (jetEta < fJetMinEta \|\| jetEta > fJetMaxEta \|\| jetPhi < fJetMinPhi \|\| jetPhi > fJetMaxPhi) {
	365	AliDebug(11,"Cut rejecting jet: Acceptance");
	366	fRejectionReason \|= kAcceptanceCut;
6c3dc2d6	367	return kFALSE;
6c3dc2d6	368	}
3c9775d9	369
68e03fc3	370	if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap) {
	371	AliDebug(11,"Cut rejecting jet: Bit map");
	372	fRejectionReason \|= kBitMapCut;
	373	return kFALSE;
	374	}
	375
6c3dc2d6	376	if (jet->Area() <= fJetAreaCut) {
6c3dc2d6	377	AliDebug(11,"Cut rejecting jet: Area");
6504c17f	378	fRejectionReason \|= kAreaCut;
6c3dc2d6	379	return kFALSE;
6c3dc2d6	380	}
3c9775d9	381
6c3dc2d6	382	if (jet->AreaEmc() < fAreaEmcCut) {
6c3dc2d6	383	AliDebug(11,"Cut rejecting jet: AreaEmc");
6504c17f	384	fRejectionReason \|= kAreaEmcCut;
6c3dc2d6	385	return kFALSE;
6c3dc2d6	386	}
3c9775d9	387
6c3dc2d6	388	if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut) {
6c3dc2d6	389	AliDebug(11,"Cut rejecting jet: ZLeading");
6504c17f	390	fRejectionReason \|= kZLeadingChCut;
6c3dc2d6	391	return kFALSE;
6c3dc2d6	392	}
3c9775d9	393
6c3dc2d6	394	if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut) {
6c3dc2d6	395	AliDebug(11,"Cut rejecting jet: ZLeadEmc");
6504c17f	396	fRejectionReason \|= kZLeadingEmcCut;
6c3dc2d6	397	return kFALSE;
6c3dc2d6	398	}
3c9775d9	399
6c3dc2d6	400	if (jet->NEF() < fNEFMinCut \|\| jet->NEF() > fNEFMaxCut) {
6c3dc2d6	401	AliDebug(11,"Cut rejecting jet: NEF");
6504c17f	402	fRejectionReason \|= kNEFCut;
6c3dc2d6	403	return kFALSE;
6c3dc2d6	404	}
3c9775d9	405
6c3dc2d6	406	if (!AcceptBiasJet(jet)) {
6c3dc2d6	407	AliDebug(11,"Cut rejecting jet: Bias");
6504c17f	408	fRejectionReason \|= kMinLeadPtCut;
6c3dc2d6	409	return kFALSE;
	410	}
	411
6504c17f	412	if (jet->MaxTrackPt() > fMaxTrackPt) {
	413	AliDebug(11,"Cut rejecting jet: MaxTrackPt");
	414	fRejectionReason \|= kMaxTrackPtCut;
	415	return kFALSE;
	416
	417	}
	418
	419	if (jet->MaxClusterPt() > fMaxClusterPt) {
	420	AliDebug(11,"Cut rejecting jet: MaxClusPt");
	421	fRejectionReason \|= kMaxClusterPtCut;
6c3dc2d6	422	return kFALSE;
6c3dc2d6	423	}
3c9775d9	424
6c3dc2d6	425	if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection)) {
6c3dc2d6	426	AliDebug(11,"Cut rejecting jet: Flavour");
6504c17f	427	fRejectionReason \|= kFlavourCut;
6c3dc2d6	428	return kFALSE;
6c3dc2d6	429	}
3c9775d9	430
6c3dc2d6	431	if(fTagStatus>-1 && jet->GetTagStatus()!=fTagStatus) {
6c3dc2d6	432	AliDebug(11,"Cut rejecting jet: tag status");
6504c17f	433	fRejectionReason \|= kTagStatus;
6c3dc2d6	434	return kFALSE;
6c3dc2d6	435	}
6504c17f	436
6504c17f	437	return kTRUE;
421b12a3	438	}
	439
	440	//________________________________________________________________________
983c25fc	441	Double_t AliJetContainer::GetLeadingHadronPt(const AliEmcalJet *jet) const
421b12a3	442	{
	443	if (fLeadingHadronType == 0) // charged leading hadron
	444	return jet->MaxTrackPt();
	445	else if (fLeadingHadronType == 1) // neutral leading hadron
	446	return jet->MaxClusterPt();
	447	else // charged or neutral
	448	return jet->MaxPartPt();
	449	}
	450
6421eeb0	451	//________________________________________________________________________
983c25fc	452	void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
6421eeb0	453	{
	454	Double_t maxClusterPt = 0;
	455	Double_t maxClusterEta = 0;
	456	Double_t maxClusterPhi = 0;
	457
	458	Double_t maxTrackPt = 0;
	459	Double_t maxTrackEta = 0;
	460	Double_t maxTrackPhi = 0;
	461
	462	if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 \|\| fLeadingHadronType == 2)) {
	463	AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
	464	if (cluster) {
	465	TLorentzVector nPart;
	466	cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
	467
	468	maxClusterEta = nPart.Eta();
	469	maxClusterPhi = nPart.Phi();
	470	maxClusterPt = nPart.Pt();
	471	}
	472	}
	473
	474	if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 \|\| fLeadingHadronType == 2)) {
	475	AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
	476	if (track) {
	477	maxTrackEta = track->Eta();
	478	maxTrackPhi = track->Phi();
	479	maxTrackPt = track->Pt();
	480	}
	481	}
	482
	483	if (maxTrackPt > maxClusterPt)
	484	mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
	485	else
	486	mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
	487	}
421b12a3	488
413ab3a3	489	//________________________________________________________________________
983c25fc	490	Double_t AliJetContainer::GetZLeadingEmc(const AliEmcalJet *jet) const
413ab3a3	491	{
413ab3a3	492
ad5d2c11	493	if (fClusterContainer && fClusterContainer->GetArray()) {
413ab3a3	494	TLorentzVector mom;
	495
	496	AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
	497	if (cluster) {
	498	cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
	499
	500	return GetZ(jet,mom);
	501	}
	502	else
	503	return -1;
	504	}
	505	else
	506	return -1;
	507	}
	508
	509	//________________________________________________________________________
983c25fc	510	Double_t AliJetContainer::GetZLeadingCharged(const AliEmcalJet *jet) const
413ab3a3	511	{
	512
	513	if (fParticleContainer && fParticleContainer->GetArray() ) {
	514	TLorentzVector mom;
	515
	516	AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
	517	if (track) {
	518	mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
	519
	520	return GetZ(jet,mom);
	521	}
	522	else
	523	return -1;
	524	}
	525	else
	526	return -1;
	527	}
	528
	529	//________________________________________________________________________
983c25fc	530	Double_t AliJetContainer::GetZ(const AliEmcalJet *jet, TLorentzVector mom) const
413ab3a3	531	{
	532
	533	Double_t pJetSq = jet->Px()jet->Px() + jet->Py()jet->Py() + jet->Pz()*jet->Pz();
	534
ad5d2c11	535	if(pJetSq>1e-6)
413ab3a3	536	return (mom.Px()jet->Px() + mom.Py()jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
	537	else {
	538	AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
	539	return -1;
	540	}
	541
	542	}
	543
421b12a3	544	//________________________________________________________________________
	545	void AliJetContainer::SetJetEtaPhiEMCAL()
	546	{
	547	//Set default cuts for full jets
	548
	549	if(!fGeom) SetEMCALGeometry();
	550	if(fGeom) {
	551	SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
56f370b0	552
56f370b0	553	if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
9dbcdad7	554	SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
56f370b0	555	else
56f370b0	556	SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
421b12a3	557	}
421b12a3	558	else {
56f370b0	559	AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
421b12a3	560	SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
254370e1	561	SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
421b12a3	562	}
421b12a3	563	}
	564
	565	//________________________________________________________________________
	566	void AliJetContainer::SetJetEtaPhiTPC()
	567	{
56f370b0	568	//Set default cuts for charged jets
421b12a3	569
56f370b0	570	SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
421b12a3	571	SetJetPhiLimits(-10, 10);
421b12a3	572	}
421b12a3	573
6c3dc2d6	574	//________________________________________________________________________
	575	void AliJetContainer::PrintCuts()
	576	{
	577	// Reset cuts to default values
	578
	579	Printf("PtBiasJetTrack: %f",fPtBiasJetTrack);
	580	Printf("PtBiasJetClus: %f",fPtBiasJetClus);
	581	Printf("JetPtCut: %f", fJetPtCut);
	582	Printf("JetAreaCut: %f",fJetAreaCut);
	583	Printf("AreaEmcCut: %f",fAreaEmcCut);
	584	Printf("JetMinEta: %f", fJetMinEta);
	585	Printf("JetMaxEta: %f", fJetMaxEta);
	586	Printf("JetMinPhi: %f", fJetMinPhi);
	587	Printf("JetMaxPhi: %f", fJetMaxPhi);
	588	Printf("MaxClusterPt: %f",fMaxClusterPt);
	589	Printf("MaxTrackPt: %f",fMaxTrackPt);
	590	Printf("LeadingHadronType: %d",fLeadingHadronType);
	591	Printf("ZLeadingEmcCut: %f",fZLeadingEmcCut);
	592	Printf("ZLeadingChCut: %f",fZLeadingChCut);
	593
	594	}
	595
421b12a3	596	//________________________________________________________________________
	597	void AliJetContainer::ResetCuts()
	598	{
	599	// Reset cuts to default values
	600
	601	fPtBiasJetTrack = 0;
	602	fPtBiasJetClus = 0;
	603	fJetPtCut = 1;
	604	fJetAreaCut = -1;
	605	fAreaEmcCut = 0;
	606	fJetMinEta = -0.9;
	607	fJetMaxEta = 0.9;
	608	fJetMinPhi = -10;
	609	fJetMaxPhi = 10;
	610	fMaxClusterPt = 1000;
	611	fMaxTrackPt = 100;
	612	fLeadingHadronType = 0;
413ab3a3	613	fZLeadingEmcCut = 10.;
413ab3a3	614	fZLeadingChCut = 10.;
421b12a3	615	}
b6f970ad	616
	617	//________________________________________________________________________
	618	void AliJetContainer::SetClassName(const char *clname)
	619	{
	620	// Set the class name
	621
	622	TClass cls(clname);
	623	if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
	624	else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
	625	}
ba7663ad	626
ba7663ad	627	//________________________________________________________________________
1705bab1	628	Double_t AliJetContainer::GetFractionSharedPt(const AliEmcalJet *jet1) const
ba7663ad	629	{
ba7663ad	630	//
f081fbfd	631	// Get fraction of shared pT between matched jets
f081fbfd	632	// Uses ClosestJet() jet pT as baseline: fraction = \Sum_{const,jet1} pT,const,i / pT,jet,closest
ba7663ad	633	// Only works if tracks array of both jets is the same
	634	//
	635
	636	AliEmcalJet *jet2 = jet1->ClosestJet();
	637	if(!jet2) return -1;
	638
	639	Double_t fraction = 0.;
	640	Double_t jetPt2 = jet2->Pt();
	641
	642	if(jetPt2>0) {
	643	Double_t sumPt = 0.;
	644	AliVParticle *vpf = 0x0;
	645	Int_t iFound = 0;
	646	for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
	647	Int_t idx = (Int_t)jet2->TrackAt(icc);
	648	iFound = 0;
	649	for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
	650	if(idx == jet1->TrackAt(icf) && iFound==0 ) {
	651	iFound=1;
	652	vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fParticleContainer->GetArray()));
	653	if(vpf) sumPt += vpf->Pt();
	654	continue;
	655	}
	656	}
	657	}
	658	fraction = sumPt/jetPt2;
	659	} else
	660	fraction = -1;
	661
	662	return fraction;
	663	}
	664