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

//
// Emcal jet class.
//
// Author: C.Loizides

#include "AliEmcalJet.h"

#include "AliLog.h"
#include "Riostream.h"

ClassImp(AliEmcalJet)

//__________________________________________________________________________________________________
AliEmcalJet::AliEmcalJet() : 
  AliVParticle(), 
  fPt(0), 
  fEta(0), 
  fPhi(0), 
  fM(0), 
  fNEF(0),
  fArea(0),       
  fAreaEta(0),       
  fAreaPhi(0),       
  fAreaEmc(-1), 
  fAxisInEmcal(0), 
  fFlavourTagging(0),
  fMaxCPt(0), 
  fMaxNPt(0), 
  fMCPt(0),
  fNn(0), 
  fNch(0),        
  fPtEmc(0),
  fNEmc(0),
  fClusterIDs(),
  fTrackIDs(),
  fMatched(2),
  fMatchingType(0),
  fTaggedJet(0x0),
  fTagStatus(-1),
  fPtSub(0),
  fPtVectSub(0),
  fTriggers(0),
  fJetShapeMassFirstDer(0),
  fJetShapeMassSecondDer(0),
  fJetShapeMassFirstSub(0),
  fJetShapeMassSecondSub(0),
  fLabel(-1),
  fGRNumerator(0),
  fGRDenominator(0),
  fGRNumeratorSub(0),
  fGRDenominatorSub(0),
  fJetShapeAngularityFirstDer(0),
  fJetShapeAngularitySecondDer(0),
  fJetShapeAngularityFirstSub(0),
  fJetShapeAngularitySecondSub(0),
  fJetShapepTDFirstDer(0),
  fJetShapepTDSecondDer(0),
  fJetShapepTDFirstSub(0),
  fJetShapepTDSecondSub(0),
  fJetShapeCircularityFirstDer(0),
  fJetShapeCircularitySecondDer(0),
  fJetShapeCircularityFirstSub(0),
  fJetShapeCircularitySecondSub(0),
  fJetShapeConstituentFirstDer(0),
  fJetShapeConstituentSecondDer(0),
  fJetShapeConstituentFirstSub(0),
  fJetShapeConstituentSecondSub(0),
  fJetShapeLeSubFirstDer(0),
  fJetShapeLeSubSecondDer(0),
  fJetShapeLeSubFirstSub(0),
  fJetShapeLeSubSecondSub(0)
{
  // Constructor.
  fClosestJets[0] = 0;
  fClosestJets[1] = 0; 
  fClosestJetsDist[0] = 999; 
  fClosestJetsDist[1] = 999; 
}

//__________________________________________________________________________________________________
AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) : 
  AliVParticle(), 
  fPt(TMath::Sqrt(px*px+py*py)), 
  fEta(TMath::ASinH(pz/fPt)),
  fPhi(0), 
  fM(0), 
  fNEF(0), 
  fArea(0), 
  fAreaEta(0),       
  fAreaPhi(0),       
  fAreaEmc(-1), 
  fAxisInEmcal(0),
  fFlavourTagging(0),
  fMaxCPt(0), 
  fMaxNPt(0), 
  fMCPt(0),
  fNn(0),
  fNch(0),
  fPtEmc(0),
  fNEmc(0),
  fClusterIDs(), 
  fTrackIDs(),
  fMatched(2),
  fMatchingType(0),
  fTaggedJet(0x0),
  fTagStatus(-1),
  fPtSub(0),
  fPtVectSub(0),
  fTriggers(0),
  fJetShapeMassFirstDer(0),
  fJetShapeMassSecondDer(0),
  fJetShapeMassFirstSub(0),
  fJetShapeMassSecondSub(0),
  fLabel(-1),
  fGRNumerator(0),
  fGRDenominator(0),
  fGRNumeratorSub(0),
  fGRDenominatorSub(0),
  fJetShapeAngularityFirstDer(0),
  fJetShapeAngularitySecondDer(0),
  fJetShapeAngularityFirstSub(0),
  fJetShapeAngularitySecondSub(0),
  fJetShapepTDFirstDer(0),
  fJetShapepTDSecondDer(0),
  fJetShapepTDFirstSub(0),
  fJetShapepTDSecondSub(0),
  fJetShapeCircularityFirstDer(0),
  fJetShapeCircularitySecondDer(0),
  fJetShapeCircularityFirstSub(0),
  fJetShapeCircularitySecondSub(0),
  fJetShapeConstituentFirstDer(0),
  fJetShapeConstituentSecondDer(0),
  fJetShapeConstituentFirstSub(0),
  fJetShapeConstituentSecondSub(0),
   fJetShapeLeSubFirstDer(0),
  fJetShapeLeSubSecondDer(0),
  fJetShapeLeSubFirstSub(0),
  fJetShapeLeSubSecondSub(0)
{    
  // Constructor.

  if (fPt != 0) {
    fPhi = TMath::ATan2(py, px);
    if (fPhi<0.) 
      fPhi += 2. * TMath::Pi();
  }

  fClosestJets[0] = 0; 
  fClosestJets[1] = 0;
  fClosestJetsDist[0] = 999; 
  fClosestJetsDist[1] = 999;
}

//_________________________________________________________________________________________________
AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
  AliVParticle(), 
  fPt(pt), 
  fEta(eta), 
  fPhi(phi), 
  fM(m), 
  fNEF(0), 
  fArea(0), 
  fAreaEta(0),       
  fAreaPhi(0),       
  fAreaEmc(-1), 
  fAxisInEmcal(0),
  fFlavourTagging(0),
  fMaxCPt(0), 
  fMaxNPt(0),
  fMCPt(0),
  fNn(0),
  fNch(0), 
  fPtEmc(0),
  fNEmc(0),
  fClusterIDs(), 
  fTrackIDs(),
  fMatched(2),
  fMatchingType(0),
  fTaggedJet(0x0),
  fTagStatus(-1),
  fPtSub(0),
  fPtVectSub(0),
  fTriggers(0),
  fJetShapeMassFirstDer(0),
  fJetShapeMassSecondDer(0),
  fJetShapeMassFirstSub(0),
  fJetShapeMassSecondSub(0),
  fLabel(-1),
  fGRNumerator(0),
  fGRDenominator(0),
  fGRNumeratorSub(0),
  fGRDenominatorSub(0),
  fJetShapeAngularityFirstDer(0),
  fJetShapeAngularitySecondDer(0),
  fJetShapeAngularityFirstSub(0),
  fJetShapeAngularitySecondSub(0),
  fJetShapepTDFirstDer(0),
  fJetShapepTDSecondDer(0),
  fJetShapepTDFirstSub(0),
  fJetShapepTDSecondSub(0),
  fJetShapeCircularityFirstDer(0),
  fJetShapeCircularitySecondDer(0),
  fJetShapeCircularityFirstSub(0),
  fJetShapeCircularitySecondSub(0),
  fJetShapeConstituentFirstDer(0),
  fJetShapeConstituentSecondDer(0),
  fJetShapeConstituentFirstSub(0),
  fJetShapeConstituentSecondSub(0),
  fJetShapeLeSubFirstDer(0),
  fJetShapeLeSubSecondDer(0),
  fJetShapeLeSubFirstSub(0),
  fJetShapeLeSubSecondSub(0)

{
  // Constructor.

  if (fPhi<0.) 
    fPhi += TMath::TwoPi();

  fClosestJets[0] = 0; 
  fClosestJets[1] = 0;
  fClosestJetsDist[0] = 999; 
  fClosestJetsDist[1] = 999;
}

//_________________________________________________________________________________________________
AliEmcalJet::AliEmcalJet(const AliEmcalJet &jet) :
  AliVParticle(jet),
  fPt(jet.fPt), 
  fEta(jet.fEta), 
  fPhi(jet.fPhi), 
  fM(jet.fM), 
  fNEF(jet.fNEF), 
  fArea(jet.fArea), 
  fAreaEta(jet.fAreaEta),       
  fAreaPhi(jet.fAreaPhi),       
  fAreaEmc(jet.fAreaEmc), 
  fAxisInEmcal(jet.fAxisInEmcal),
  fFlavourTagging(jet.fFlavourTagging),
  fMaxCPt(jet.fMaxCPt), 
  fMaxNPt(jet.fMaxNPt), 
  fMCPt(jet.fMCPt),
  fNn(jet.fNn),
  fNch(jet.fNch),
  fPtEmc(jet.fPtEmc),
  fNEmc(jet.fNEmc),
  fClusterIDs(jet.fClusterIDs), 
  fTrackIDs(jet.fTrackIDs),
  fMatched(jet.fMatched),
  fMatchingType(jet.fMatchingType),
  fTaggedJet(jet.fTaggedJet),
  fTagStatus(jet.fTagStatus),
  fPtSub(jet.fPtSub),
  fPtVectSub(jet.fPtVectSub),
  fTriggers(jet.fTriggers),
  fJetShapeMassFirstDer(jet.fJetShapeMassFirstDer),
  fJetShapeMassSecondDer(jet.fJetShapeMassSecondDer),
  fJetShapeMassFirstSub(jet.fJetShapeMassFirstSub),
  fJetShapeMassSecondSub(jet.fJetShapeMassSecondSub),
  fLabel(jet.fLabel),
  fGRNumerator(jet.fGRNumerator),
  fGRDenominator(jet.fGRDenominator),
  fGRNumeratorSub(jet.fGRNumeratorSub),
  fGRDenominatorSub(jet.fGRDenominatorSub),
  fJetShapeAngularityFirstDer(jet.fJetShapeAngularityFirstDer),
  fJetShapeAngularitySecondDer(jet.fJetShapeAngularitySecondDer),
  fJetShapeAngularityFirstSub(jet.fJetShapeAngularityFirstSub),
  fJetShapeAngularitySecondSub(jet.fJetShapeAngularitySecondSub),
  fJetShapepTDFirstDer(jet.fJetShapepTDFirstDer),
  fJetShapepTDSecondDer(jet.fJetShapepTDSecondDer),
  fJetShapepTDFirstSub(jet.fJetShapepTDFirstSub),
  fJetShapepTDSecondSub(jet.fJetShapepTDSecondSub),
  fJetShapeCircularityFirstDer(jet.fJetShapeCircularityFirstDer),
  fJetShapeCircularitySecondDer(jet.fJetShapeCircularitySecondDer),
  fJetShapeCircularityFirstSub(jet.fJetShapeCircularityFirstSub),
  fJetShapeCircularitySecondSub(jet.fJetShapeCircularitySecondSub),
  fJetShapeConstituentFirstDer(jet.fJetShapeConstituentFirstDer),
  fJetShapeConstituentSecondDer(jet.fJetShapeConstituentSecondDer),
  fJetShapeConstituentFirstSub(jet.fJetShapeConstituentFirstSub),
  fJetShapeConstituentSecondSub(jet.fJetShapeConstituentSecondSub),
  fJetShapeLeSubFirstDer(jet.fJetShapeLeSubFirstDer),
  fJetShapeLeSubSecondDer(jet.fJetShapeLeSubSecondDer),
  fJetShapeLeSubFirstSub(jet.fJetShapeLeSubFirstSub),
  fJetShapeLeSubSecondSub(jet.fJetShapeLeSubSecondSub)
{
  // Copy constructor.
  fClosestJets[0]     = jet.fClosestJets[0]; 
  fClosestJets[1]     = jet.fClosestJets[1]; 
  fClosestJetsDist[0] = jet.fClosestJetsDist[0];  
  fClosestJetsDist[1] = jet.fClosestJetsDist[1]; 
}

//_________________________________________________________________________________________________
AliEmcalJet &AliEmcalJet::operator=(const AliEmcalJet &jet)
{
  // Assignment operator.

  if (this!=&jet) {
    AliVParticle::operator=(jet);
    fPt                 = jet.fPt;
    fEta                = jet.fEta;
    fPhi                = jet.fPhi;
    fM                  = jet.fM; 
    fNEF                = jet.fNEF;
    fArea               = jet.fArea; 
    fAreaEta            = jet.fAreaEta; 
    fAreaPhi            = jet.fAreaPhi; 
    fAreaEmc            = jet.fAreaEmc; 
    fAxisInEmcal        = jet.fAxisInEmcal; 
    fFlavourTagging     = jet.fFlavourTagging;
    fMaxCPt             = jet.fMaxCPt; 
    fMaxNPt             = jet.fMaxNPt;
    fMCPt               = jet.fMCPt;
    fNn                 = jet.fNn;
    fNch                = jet.fNch;
    fPtEmc              = jet.fPtEmc;
    fNEmc               = jet.fNEmc;
    fClusterIDs         = jet.fClusterIDs;
    fTrackIDs           = jet.fTrackIDs;
    fClosestJets[0]     = jet.fClosestJets[0]; 
    fClosestJets[1]     = jet.fClosestJets[1]; 
    fClosestJetsDist[0] = jet.fClosestJetsDist[0];  
    fClosestJetsDist[1] = jet.fClosestJetsDist[1]; 
    fMatched            = jet.fMatched;
    fTaggedJet          = jet.fTaggedJet;
    fTagStatus          = jet.fTagStatus;
    fPtSub              = jet.fPtSub;
    fPtVectSub          = jet.fPtVectSub;
    fTriggers           = jet.fTriggers;
    fJetShapeMassFirstDer  = jet.fJetShapeMassFirstDer;
    fJetShapeMassSecondDer = jet.fJetShapeMassSecondDer;
    fJetShapeMassFirstSub  = jet.fJetShapeMassFirstSub;
    fJetShapeMassSecondSub = jet.fJetShapeMassSecondSub;
    fLabel              = jet.fLabel;
    fGRNumerator        = jet.fGRNumerator;
    fGRDenominator      = jet.fGRDenominator;
    fGRNumeratorSub     = jet.fGRNumeratorSub;
    fGRDenominatorSub   = jet.fGRDenominatorSub;
    fJetShapeAngularityFirstDer  = jet.fJetShapeAngularityFirstDer;
    fJetShapeAngularitySecondDer = jet.fJetShapeAngularitySecondDer;
    fJetShapeAngularityFirstSub  = jet.fJetShapeAngularityFirstSub;
    fJetShapeAngularitySecondSub = jet.fJetShapeAngularitySecondSub;
    fJetShapepTDFirstDer  = jet.fJetShapepTDFirstDer;
    fJetShapepTDSecondDer = jet.fJetShapepTDSecondDer;
    fJetShapepTDFirstSub  = jet.fJetShapepTDFirstSub;
    fJetShapepTDSecondSub = jet.fJetShapepTDSecondSub;
    fJetShapeCircularityFirstDer  = jet.fJetShapeCircularityFirstDer;
    fJetShapeCircularitySecondDer = jet.fJetShapeCircularitySecondDer;
    fJetShapeCircularityFirstSub  = jet.fJetShapeCircularityFirstSub;
    fJetShapeCircularitySecondSub = jet.fJetShapeCircularitySecondSub;
    fJetShapeConstituentFirstDer  = jet.fJetShapeConstituentFirstDer;
    fJetShapeConstituentSecondDer = jet.fJetShapeConstituentSecondDer;
    fJetShapeConstituentFirstSub  = jet.fJetShapeConstituentFirstSub;
    fJetShapeConstituentSecondSub = jet.fJetShapeConstituentSecondSub;
    fJetShapeLeSubFirstDer  = jet.fJetShapeLeSubFirstDer;
    fJetShapeLeSubSecondDer = jet.fJetShapeLeSubSecondDer;
    fJetShapeLeSubFirstSub  = jet.fJetShapeLeSubFirstSub;
    fJetShapeLeSubSecondSub = jet.fJetShapeLeSubSecondSub;
  }

  return *this;
}

//_________________________________________________________________________________________________
Int_t AliEmcalJet::Compare(const TObject* obj) const
{
  //Return -1 if this is smaller than obj, 0 if objects are equal and 1 if this is larger than obj.

  const AliEmcalJet *jet = static_cast<const AliEmcalJet *>(obj);
  if (!obj)
    return 0;
  if (Pt()>jet->Pt())
    return -1;
  return 1;
}

//__________________________________________________________________________________________________
void AliEmcalJet::GetMom(TLorentzVector &vec) const
{
  // Return momentum as four vector.

  Double_t p = fPt *TMath::CosH(fEta);
  vec.SetPtEtaPhiE(fPt,fEta,fPhi,TMath::Sqrt(p*p+fM*fM));
}

//__________________________________________________________________________________________________
void AliEmcalJet::Print(Option_t* /*option*/) const
{
  // Print jet information.

  Printf("Jet pt=%.2f, eta=%.2f, phi=%.2f, area=%.2f, NEF=%.2f", fPt, fEta, fPhi, fArea, fNEF);
}

//__________________________________________________________________________________________________
Double_t AliEmcalJet::PtSubVect(Double_t rho) const
{
  // Return vectorial subtracted transverse momentum.

  Double_t dx = Px() - rho * fArea * TMath::Cos(fAreaPhi);
  Double_t dy = Py() - rho * fArea * TMath::Sin(fAreaPhi);
  //Double_t dz = Pz() - rho * fArea * TMath::SinH(fAreaEta);
  return TMath::Sqrt(dx*dx+dy*dy);
}

//__________________________________________________________________________________________________
void AliEmcalJet::SortConstituents()
{
  // Sort constituent by index (increasing).

  std::sort(fClusterIDs.GetArray(), fClusterIDs.GetArray() + fClusterIDs.GetSize());
  std::sort(fTrackIDs.GetArray(), fTrackIDs.GetArray() + fTrackIDs.GetSize());
}

//__________________________________________________________________________________________________
Double_t AliEmcalJet::DeltaR (const AliVParticle* part) const
{ // Helper function to calculate the distance between two jets or a jet and a particle
    Double_t dPhi = this->Phi() - part->Phi();
    Double_t dEta = this->Eta() - part->Eta();
    dPhi = TVector2::Phi_mpi_pi ( dPhi );

    return TMath::Sqrt ( dPhi * dPhi + dEta * dEta );
}


//__________________________________________________________________________________________________
std::vector<int> AliEmcalJet::SortConstituentsPt( TClonesArray *tracks ) const
{   //___________________________________________
    // Sorting by p_T (decreasing) jet constituents
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    typedef std::pair<Double_t, Int_t> ptidx_pair;

    // Create vector for Pt sorting
    std::vector<ptidx_pair> pair_list ;

    for ( Int_t i_entry = 0; i_entry < GetNumberOfTracks(); i_entry++ )
        {
        AliVParticle *track = TrackAt(i_entry, tracks);
        if (!track)
            {
            AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)", i_entry, tracks->GetName(), TrackAt(i_entry), tracks->GetEntriesFast()));
            continue;
            }

        pair_list.push_back( std::make_pair ( track->Pt(), i_entry ) );
        }

    std::stable_sort( pair_list.begin() , pair_list.end() , sort_descend() );

    // return an vector of indexes of constituents (sorted descending by pt)
    std::vector <int> index_sorted_list;

    for ( std::vector< std::pair<Double_t,Int_t> >::iterator it = pair_list.begin(); it != pair_list.end(); ++it)
        { index_sorted_list.push_back( (*it).second ); } // populating the return object with indexes of sorted tracks

    return index_sorted_list;
}

//________________________________________________________________________
Double_t AliEmcalJet::GetZ ( const Double_t trkPx, const Double_t trkPy, const Double_t trkPz ) const
    {
    // Get the z of a constituent inside of a jet
    Double_t pJetSq = this->Px() * this->Px() + this->Py() * this->Py() + this->Pz() * this->Pz();

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

    }

//________________________________________________________________________
Double_t AliEmcalJet::GetZ ( const AliVParticle* trk )          const
    {
    // Get Z of constituent trk
    return GetZ ( trk->Px(), trk->Py(), trk->Pz() );
    }

//__________________________________________________________________________________________________
AliVParticle* AliEmcalJet::GetLeadingTrack(TClonesArray *tracks) const
{
  AliVParticle* maxTrack = 0;
  for (Int_t i = 0; i < GetNumberOfTracks(); i++) {
    AliVParticle *track = TrackAt(i, tracks);
    if (!track) {
      AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)",
		    i,tracks->GetName(),TrackAt(i),tracks->GetEntriesFast()));
      continue;
    }
    if (!maxTrack || track->Pt() > maxTrack->Pt()) 
      maxTrack = track;
  }

  return maxTrack;
}

//__________________________________________________________________________________________________
AliVCluster* AliEmcalJet::GetLeadingCluster(TClonesArray *clusters) const
{
  AliVCluster* maxCluster = 0;
  for (Int_t i = 0; i < GetNumberOfClusters(); i++) {
    AliVCluster *cluster = ClusterAt(i, clusters);
    if (!cluster) {
      AliError(Form("Unable to find jet cluster %d in collection %s (pos in collection %d, max %d)",
		    i,clusters->GetName(),ClusterAt(i),clusters->GetEntriesFast()));
      continue;
    }
    if (!maxCluster || cluster->E() > maxCluster->E()) 
      maxCluster = cluster;
  }

  return maxCluster;
}

//__________________________________________________________________________________________________
void AliEmcalJet::ResetMatching()
{
  fClosestJets[0] = 0;
  fClosestJets[1] = 0; 
  fClosestJetsDist[0] = 999; 
  fClosestJetsDist[1] = 999; 
  fMatched = 2;
}

//__________________________________________________________________________________________________
void AliEmcalJet::PrintGR() {
  for(Int_t i = 0; i<fGRNumerator.GetSize(); i++) {
    Printf("num[%d] = %f",i,fGRNumerator.At(i));
  }
}
Commit	Line	Data
04f6e9de	1	//
	2	// Emcal jet class.
	3	//
cd231d42	4	// Author: C.Loizides
7efbea04	5
914d486c	6	#include "AliEmcalJet.h"
7df864a3	7
59ebe296	8	#include "AliLog.h"
	9	#include "Riostream.h"
	10
914d486c	11	ClassImp(AliEmcalJet)
7efbea04	12
96919f91	13	//__________________________________________________________________________________________________
	14	AliEmcalJet::AliEmcalJet() :
	15	AliVParticle(),
	16	fPt(0),
	17	fEta(0),
	18	fPhi(0),
	19	fM(0),
	20	fNEF(0),
	21	fArea(0),
db8a0184	22	fAreaEta(0),
db8a0184	23	fAreaPhi(0),
96919f91	24	fAreaEmc(-1),
96919f91	25	fAxisInEmcal(0),
3c9775d9	26	fFlavourTagging(0),
96919f91	27	fMaxCPt(0),
	28	fMaxNPt(0),
	29	fMCPt(0),
	30	fNn(0),
	31	fNch(0),
8895e61b	32	fPtEmc(0),
8895e61b	33	fNEmc(0),
96919f91	34	fClusterIDs(),
96919f91	35	fTrackIDs(),
8895e61b	36	fMatched(2),
cfc2ac24	37	fMatchingType(0),
254370e1	38	fTaggedJet(0x0),
254370e1	39	fTagStatus(-1),
db8a0184	40	fPtSub(0),
85d48773	41	fPtVectSub(0),
8082a80b	42	fTriggers(0),
	43	fJetShapeMassFirstDer(0),
	44	fJetShapeMassSecondDer(0),
	45	fJetShapeMassFirstSub(0),
	46	fJetShapeMassSecondSub(0),
726c9488	47	fLabel(-1),
	48	fGRNumerator(0),
	49	fGRDenominator(0),
	50	fGRNumeratorSub(0),
0d13a63c	51	fGRDenominatorSub(0),
	52	fJetShapeAngularityFirstDer(0),
	53	fJetShapeAngularitySecondDer(0),
	54	fJetShapeAngularityFirstSub(0),
	55	fJetShapeAngularitySecondSub(0),
	56	fJetShapepTDFirstDer(0),
	57	fJetShapepTDSecondDer(0),
	58	fJetShapepTDFirstSub(0),
	59	fJetShapepTDSecondSub(0),
	60	fJetShapeCircularityFirstDer(0),
	61	fJetShapeCircularitySecondDer(0),
	62	fJetShapeCircularityFirstSub(0),
	63	fJetShapeCircularitySecondSub(0),
	64	fJetShapeConstituentFirstDer(0),
	65	fJetShapeConstituentSecondDer(0),
	66	fJetShapeConstituentFirstSub(0),
	67	fJetShapeConstituentSecondSub(0),
	68	fJetShapeLeSubFirstDer(0),
	69	fJetShapeLeSubSecondDer(0),
	70	fJetShapeLeSubFirstSub(0),
	71	fJetShapeLeSubSecondSub(0)
96919f91	72	{
96919f91	73	// Constructor.
96919f91	74	fClosestJets[0] = 0;
	75	fClosestJets[1] = 0;
	76	fClosestJetsDist[0] = 999;
	77	fClosestJetsDist[1] = 999;
96919f91	78	}
96919f91	79
7efbea04	80	//__________________________________________________________________________________________________
101cefde	81	AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
	82	AliVParticle(),
	83	fPt(TMath::Sqrt(pxpx+pypy)),
	84	fEta(TMath::ASinH(pz/fPt)),
	85	fPhi(0),
	86	fM(0),
	87	fNEF(0),
	88	fArea(0),
db8a0184	89	fAreaEta(0),
db8a0184	90	fAreaPhi(0),
101cefde	91	fAreaEmc(-1),
101cefde	92	fAxisInEmcal(0),
3c9775d9	93	fFlavourTagging(0),
101cefde	94	fMaxCPt(0),
101cefde	95	fMaxNPt(0),
96919f91	96	fMCPt(0),
	97	fNn(0),
	98	fNch(0),
8895e61b	99	fPtEmc(0),
8895e61b	100	fNEmc(0),
101cefde	101	fClusterIDs(),
96919f91	102	fTrackIDs(),
8895e61b	103	fMatched(2),
cfc2ac24	104	fMatchingType(0),
254370e1	105	fTaggedJet(0x0),
254370e1	106	fTagStatus(-1),
db8a0184	107	fPtSub(0),
85d48773	108	fPtVectSub(0),
8082a80b	109	fTriggers(0),
	110	fJetShapeMassFirstDer(0),
	111	fJetShapeMassSecondDer(0),
	112	fJetShapeMassFirstSub(0),
	113	fJetShapeMassSecondSub(0),
726c9488	114	fLabel(-1),
	115	fGRNumerator(0),
	116	fGRDenominator(0),
	117	fGRNumeratorSub(0),
0d13a63c	118	fGRDenominatorSub(0),
3ec5da8e	119	fJetShapeAngularityFirstDer(0),
0d13a63c	120	fJetShapeAngularitySecondDer(0),
	121	fJetShapeAngularityFirstSub(0),
	122	fJetShapeAngularitySecondSub(0),
3ec5da8e	123	fJetShapepTDFirstDer(0),
0d13a63c	124	fJetShapepTDSecondDer(0),
	125	fJetShapepTDFirstSub(0),
	126	fJetShapepTDSecondSub(0),
	127	fJetShapeCircularityFirstDer(0),
	128	fJetShapeCircularitySecondDer(0),
	129	fJetShapeCircularityFirstSub(0),
	130	fJetShapeCircularitySecondSub(0),
	131	fJetShapeConstituentFirstDer(0),
	132	fJetShapeConstituentSecondDer(0),
	133	fJetShapeConstituentFirstSub(0),
	134	fJetShapeConstituentSecondSub(0),
	135	fJetShapeLeSubFirstDer(0),
	136	fJetShapeLeSubSecondDer(0),
	137	fJetShapeLeSubFirstSub(0),
	138	fJetShapeLeSubSecondSub(0)
7efbea04	139	{
	140	// Constructor.
	141
	142	if (fPt != 0) {
	143	fPhi = TMath::ATan2(py, px);
	144	if (fPhi<0.)
	145	fPhi += 2. * TMath::Pi();
	146	}
96919f91	147
	148	fClosestJets[0] = 0;
	149	fClosestJets[1] = 0;
	150	fClosestJetsDist[0] = 999;
	151	fClosestJetsDist[1] = 999;
7efbea04	152	}
7efbea04	153
9cffd2d0	154	//_________________________________________________________________________________________________
914d486c	155	AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
9cffd2d0	156	AliVParticle(),
f8087a81	157	fPt(pt),
	158	fEta(eta),
	159	fPhi(phi),
	160	fM(m),
	161	fNEF(0),
	162	fArea(0),
db8a0184	163	fAreaEta(0),
db8a0184	164	fAreaPhi(0),
101cefde	165	fAreaEmc(-1),
101cefde	166	fAxisInEmcal(0),
3c9775d9	167	fFlavourTagging(0),
f8087a81	168	fMaxCPt(0),
96919f91	169	fMaxNPt(0),
	170	fMCPt(0),
	171	fNn(0),
	172	fNch(0),
8895e61b	173	fPtEmc(0),
8895e61b	174	fNEmc(0),
f8087a81	175	fClusterIDs(),
96919f91	176	fTrackIDs(),
8895e61b	177	fMatched(2),
cfc2ac24	178	fMatchingType(0),
254370e1	179	fTaggedJet(0x0),
254370e1	180	fTagStatus(-1),
db8a0184	181	fPtSub(0),
85d48773	182	fPtVectSub(0),
8082a80b	183	fTriggers(0),
	184	fJetShapeMassFirstDer(0),
	185	fJetShapeMassSecondDer(0),
	186	fJetShapeMassFirstSub(0),
	187	fJetShapeMassSecondSub(0),
726c9488	188	fLabel(-1),
	189	fGRNumerator(0),
	190	fGRDenominator(0),
	191	fGRNumeratorSub(0),
0d13a63c	192	fGRDenominatorSub(0),
	193	fJetShapeAngularityFirstDer(0),
	194	fJetShapeAngularitySecondDer(0),
	195	fJetShapeAngularityFirstSub(0),
	196	fJetShapeAngularitySecondSub(0),
	197	fJetShapepTDFirstDer(0),
	198	fJetShapepTDSecondDer(0),
	199	fJetShapepTDFirstSub(0),
	200	fJetShapepTDSecondSub(0),
	201	fJetShapeCircularityFirstDer(0),
	202	fJetShapeCircularitySecondDer(0),
	203	fJetShapeCircularityFirstSub(0),
	204	fJetShapeCircularitySecondSub(0),
	205	fJetShapeConstituentFirstDer(0),
	206	fJetShapeConstituentSecondDer(0),
	207	fJetShapeConstituentFirstSub(0),
	208	fJetShapeConstituentSecondSub(0),
	209	fJetShapeLeSubFirstDer(0),
	210	fJetShapeLeSubSecondDer(0),
	211	fJetShapeLeSubFirstSub(0),
	212	fJetShapeLeSubSecondSub(0)
	213
9cffd2d0	214	{
	215	// Constructor.
	216
4f149872	217	if (fPhi<0.)
4f149872	218	fPhi += TMath::TwoPi();
96919f91	219
	220	fClosestJets[0] = 0;
	221	fClosestJets[1] = 0;
	222	fClosestJetsDist[0] = 999;
	223	fClosestJetsDist[1] = 999;
9cffd2d0	224	}
9cffd2d0	225
7efbea04	226	//_________________________________________________________________________________________________
914d486c	227	AliEmcalJet::AliEmcalJet(const AliEmcalJet &jet) :
7efbea04	228	AliVParticle(jet),
f8087a81	229	fPt(jet.fPt),
	230	fEta(jet.fEta),
	231	fPhi(jet.fPhi),
	232	fM(jet.fM),
	233	fNEF(jet.fNEF),
	234	fArea(jet.fArea),
db8a0184	235	fAreaEta(jet.fAreaEta),
db8a0184	236	fAreaPhi(jet.fAreaPhi),
101cefde	237	fAreaEmc(jet.fAreaEmc),
101cefde	238	fAxisInEmcal(jet.fAxisInEmcal),
3c9775d9	239	fFlavourTagging(jet.fFlavourTagging),
f8087a81	240	fMaxCPt(jet.fMaxCPt),
f8087a81	241	fMaxNPt(jet.fMaxNPt),
96919f91	242	fMCPt(jet.fMCPt),
	243	fNn(jet.fNn),
	244	fNch(jet.fNch),
8895e61b	245	fPtEmc(jet.fPtEmc),
8895e61b	246	fNEmc(jet.fNEmc),
f8087a81	247	fClusterIDs(jet.fClusterIDs),
4f149872	248	fTrackIDs(jet.fTrackIDs),
8895e61b	249	fMatched(jet.fMatched),
cfc2ac24	250	fMatchingType(jet.fMatchingType),
254370e1	251	fTaggedJet(jet.fTaggedJet),
254370e1	252	fTagStatus(jet.fTagStatus),
db8a0184	253	fPtSub(jet.fPtSub),
85d48773	254	fPtVectSub(jet.fPtVectSub),
8082a80b	255	fTriggers(jet.fTriggers),
	256	fJetShapeMassFirstDer(jet.fJetShapeMassFirstDer),
	257	fJetShapeMassSecondDer(jet.fJetShapeMassSecondDer),
	258	fJetShapeMassFirstSub(jet.fJetShapeMassFirstSub),
	259	fJetShapeMassSecondSub(jet.fJetShapeMassSecondSub),
726c9488	260	fLabel(jet.fLabel),
	261	fGRNumerator(jet.fGRNumerator),
	262	fGRDenominator(jet.fGRDenominator),
	263	fGRNumeratorSub(jet.fGRNumeratorSub),
0d13a63c	264	fGRDenominatorSub(jet.fGRDenominatorSub),
3ec5da8e	265	fJetShapeAngularityFirstDer(jet.fJetShapeAngularityFirstDer),
0d13a63c	266	fJetShapeAngularitySecondDer(jet.fJetShapeAngularitySecondDer),
	267	fJetShapeAngularityFirstSub(jet.fJetShapeAngularityFirstSub),
	268	fJetShapeAngularitySecondSub(jet.fJetShapeAngularitySecondSub),
0d13a63c	269	fJetShapepTDFirstDer(jet.fJetShapepTDFirstDer),
	270	fJetShapepTDSecondDer(jet.fJetShapepTDSecondDer),
	271	fJetShapepTDFirstSub(jet.fJetShapepTDFirstSub),
	272	fJetShapepTDSecondSub(jet.fJetShapepTDSecondSub),
0d13a63c	273	fJetShapeCircularityFirstDer(jet.fJetShapeCircularityFirstDer),
	274	fJetShapeCircularitySecondDer(jet.fJetShapeCircularitySecondDer),
	275	fJetShapeCircularityFirstSub(jet.fJetShapeCircularityFirstSub),
	276	fJetShapeCircularitySecondSub(jet.fJetShapeCircularitySecondSub),
0d13a63c	277	fJetShapeConstituentFirstDer(jet.fJetShapeConstituentFirstDer),
	278	fJetShapeConstituentSecondDer(jet.fJetShapeConstituentSecondDer),
	279	fJetShapeConstituentFirstSub(jet.fJetShapeConstituentFirstSub),
	280	fJetShapeConstituentSecondSub(jet.fJetShapeConstituentSecondSub),
0d13a63c	281	fJetShapeLeSubFirstDer(jet.fJetShapeLeSubFirstDer),
	282	fJetShapeLeSubSecondDer(jet.fJetShapeLeSubSecondDer),
	283	fJetShapeLeSubFirstSub(jet.fJetShapeLeSubFirstSub),
	284	fJetShapeLeSubSecondSub(jet.fJetShapeLeSubSecondSub)
7df864a3	285	{
101cefde	286	// Copy constructor.
4f149872	287	fClosestJets[0] = jet.fClosestJets[0];
	288	fClosestJets[1] = jet.fClosestJets[1];
	289	fClosestJetsDist[0] = jet.fClosestJetsDist[0];
	290	fClosestJetsDist[1] = jet.fClosestJetsDist[1];
7efbea04	291	}
	292
	293	//_________________________________________________________________________________________________
914d486c	294	AliEmcalJet &AliEmcalJet::operator=(const AliEmcalJet &jet)
7efbea04	295	{
	296	// Assignment operator.
	297
	298	if (this!=&jet) {
	299	AliVParticle::operator=(jet);
4f149872	300	fPt = jet.fPt;
	301	fEta = jet.fEta;
	302	fPhi = jet.fPhi;
	303	fM = jet.fM;
	304	fNEF = jet.fNEF;
	305	fArea = jet.fArea;
db8a0184	306	fAreaEta = jet.fAreaEta;
db8a0184	307	fAreaPhi = jet.fAreaPhi;
4f149872	308	fAreaEmc = jet.fAreaEmc;
4f149872	309	fAxisInEmcal = jet.fAxisInEmcal;
3c9775d9	310	fFlavourTagging = jet.fFlavourTagging;
4f149872	311	fMaxCPt = jet.fMaxCPt;
	312	fMaxNPt = jet.fMaxNPt;
	313	fMCPt = jet.fMCPt;
	314	fNn = jet.fNn;
	315	fNch = jet.fNch;
8895e61b	316	fPtEmc = jet.fPtEmc;
8895e61b	317	fNEmc = jet.fNEmc;
4f149872	318	fClusterIDs = jet.fClusterIDs;
	319	fTrackIDs = jet.fTrackIDs;
	320	fClosestJets[0] = jet.fClosestJets[0];
	321	fClosestJets[1] = jet.fClosestJets[1];
	322	fClosestJetsDist[0] = jet.fClosestJetsDist[0];
	323	fClosestJetsDist[1] = jet.fClosestJetsDist[1];
	324	fMatched = jet.fMatched;
254370e1	325	fTaggedJet = jet.fTaggedJet;
254370e1	326	fTagStatus = jet.fTagStatus;
8895e61b	327	fPtSub = jet.fPtSub;
f660c2d6	328	fPtVectSub = jet.fPtVectSub;
85d48773	329	fTriggers = jet.fTriggers;
8082a80b	330	fJetShapeMassFirstDer = jet.fJetShapeMassFirstDer;
	331	fJetShapeMassSecondDer = jet.fJetShapeMassSecondDer;
	332	fJetShapeMassFirstSub = jet.fJetShapeMassFirstSub;
	333	fJetShapeMassSecondSub = jet.fJetShapeMassSecondSub;
	334	fLabel = jet.fLabel;
726c9488	335	fGRNumerator = jet.fGRNumerator;
	336	fGRDenominator = jet.fGRDenominator;
	337	fGRNumeratorSub = jet.fGRNumeratorSub;
	338	fGRDenominatorSub = jet.fGRDenominatorSub;
3ec5da8e	339	fJetShapeAngularityFirstDer = jet.fJetShapeAngularityFirstDer;
0d13a63c	340	fJetShapeAngularitySecondDer = jet.fJetShapeAngularitySecondDer;
	341	fJetShapeAngularityFirstSub = jet.fJetShapeAngularityFirstSub;
	342	fJetShapeAngularitySecondSub = jet.fJetShapeAngularitySecondSub;
0d13a63c	343	fJetShapepTDFirstDer = jet.fJetShapepTDFirstDer;
	344	fJetShapepTDSecondDer = jet.fJetShapepTDSecondDer;
	345	fJetShapepTDFirstSub = jet.fJetShapepTDFirstSub;
	346	fJetShapepTDSecondSub = jet.fJetShapepTDSecondSub;
0d13a63c	347	fJetShapeCircularityFirstDer = jet.fJetShapeCircularityFirstDer;
	348	fJetShapeCircularitySecondDer = jet.fJetShapeCircularitySecondDer;
	349	fJetShapeCircularityFirstSub = jet.fJetShapeCircularityFirstSub;
	350	fJetShapeCircularitySecondSub = jet.fJetShapeCircularitySecondSub;
0d13a63c	351	fJetShapeConstituentFirstDer = jet.fJetShapeConstituentFirstDer;
	352	fJetShapeConstituentSecondDer = jet.fJetShapeConstituentSecondDer;
	353	fJetShapeConstituentFirstSub = jet.fJetShapeConstituentFirstSub;
	354	fJetShapeConstituentSecondSub = jet.fJetShapeConstituentSecondSub;
	355	fJetShapeLeSubFirstDer = jet.fJetShapeLeSubFirstDer;
	356	fJetShapeLeSubSecondDer = jet.fJetShapeLeSubSecondDer;
	357	fJetShapeLeSubFirstSub = jet.fJetShapeLeSubFirstSub;
	358	fJetShapeLeSubSecondSub = jet.fJetShapeLeSubSecondSub;
7efbea04	359	}
	360
	361	return *this;
7df864a3	362	}
7efbea04	363
2103dc6a	364	//_________________________________________________________________________________________________
db8a0184	365	Int_t AliEmcalJet::Compare(const TObject* obj) const
	366	{
	367	//Return -1 if this is smaller than obj, 0 if objects are equal and 1 if this is larger than obj.
	368
	369	const AliEmcalJet jet = static_cast<const AliEmcalJet >(obj);
	370	if (!obj)
	371	return 0;
	372	if (Pt()>jet->Pt())
	373	return -1;
	374	return 1;
	375	}
	376
7efbea04	377	//__________________________________________________________________________________________________
914d486c	378	void AliEmcalJet::GetMom(TLorentzVector &vec) const
7efbea04	379	{
	380	// Return momentum as four vector.
	381
	382	Double_t p = fPt *TMath::CosH(fEta);
	383	vec.SetPtEtaPhiE(fPt,fEta,fPhi,TMath::Sqrt(pp+fMfM));
	384	}
	385
	386	//__________________________________________________________________________________________________
914d486c	387	void AliEmcalJet::Print(Option_t* /option/) const
7efbea04	388	{
	389	// Print jet information.
	390
59ebe296	391	Printf("Jet pt=%.2f, eta=%.2f, phi=%.2f, area=%.2f, NEF=%.2f", fPt, fEta, fPhi, fArea, fNEF);
7efbea04	392	}
629e03fc	393
db8a0184	394	//__________________________________________________________________________________________________
	395	Double_t AliEmcalJet::PtSubVect(Double_t rho) const
	396	{
	397	// Return vectorial subtracted transverse momentum.
	398
	399	Double_t dx = Px() - rho * fArea * TMath::Cos(fAreaPhi);
	400	Double_t dy = Py() - rho * fArea * TMath::Sin(fAreaPhi);
	401	//Double_t dz = Pz() - rho * fArea * TMath::SinH(fAreaEta);
	402	return TMath::Sqrt(dxdx+dydy);
	403	}
	404
629e03fc	405	//__________________________________________________________________________________________________
	406	void AliEmcalJet::SortConstituents()
	407	{
f8087a81	408	// Sort constituent by index (increasing).
	409
	410	std::sort(fClusterIDs.GetArray(), fClusterIDs.GetArray() + fClusterIDs.GetSize());
	411	std::sort(fTrackIDs.GetArray(), fTrackIDs.GetArray() + fTrackIDs.GetSize());
629e03fc	412	}
f660c2d6	413
d1f0108e	414	//__________________________________________________________________________________________________
	415	Double_t AliEmcalJet::DeltaR (const AliVParticle* part) const
	416	{ // Helper function to calculate the distance between two jets or a jet and a particle
	417	Double_t dPhi = this->Phi() - part->Phi();
	418	Double_t dEta = this->Eta() - part->Eta();
	419	dPhi = TVector2::Phi_mpi_pi ( dPhi );
	420
	421	return TMath::Sqrt ( dPhi * dPhi + dEta * dEta );
	422	}
	423
	424
	425	//__________________________________________________________________________________________________
	426	std::vector<int> AliEmcalJet::SortConstituentsPt( TClonesArray *tracks ) const
	427	{ //___________________________________________
	428	// Sorting by p_T (decreasing) jet constituents
	429	//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	430	typedef std::pair<Double_t, Int_t> ptidx_pair;
	431
	432	// Create vector for Pt sorting
	433	std::vector<ptidx_pair> pair_list ;
	434
	435	for ( Int_t i_entry = 0; i_entry < GetNumberOfTracks(); i_entry++ )
	436	{
	437	AliVParticle *track = TrackAt(i_entry, tracks);
	438	if (!track)
	439	{
	440	AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)", i_entry, tracks->GetName(), TrackAt(i_entry), tracks->GetEntriesFast()));
	441	continue;
	442	}
	443
	444	pair_list.push_back( std::make_pair ( track->Pt(), i_entry ) );
	445	}
	446
	447	std::stable_sort( pair_list.begin() , pair_list.end() , sort_descend() );
	448
	449	// return an vector of indexes of constituents (sorted descending by pt)
	450	std::vector <int> index_sorted_list;
	451
	452	for ( std::vector< std::pair<Double_t,Int_t> >::iterator it = pair_list.begin(); it != pair_list.end(); ++it)
	453	{ index_sorted_list.push_back( (*it).second ); } // populating the return object with indexes of sorted tracks
	454
	455	return index_sorted_list;
	456	}
	457
2c932584	458	//________________________________________________________________________
	459	Double_t AliEmcalJet::GetZ ( const Double_t trkPx, const Double_t trkPy, const Double_t trkPz ) const
	460	{
	461	// Get the z of a constituent inside of a jet
	462	Double_t pJetSq = this->Px() * this->Px() + this->Py() * this->Py() + this->Pz() * this->Pz();
	463
	464	if(pJetSq>1e-6)
	465	{ return ( trkPx * this->Px() + trkPy * this->Py() + trkPz * this->Pz() ) / pJetSq ; }
	466	else
	467	{ AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq)); return -1; }
	468
	469	}
	470
	471	//________________________________________________________________________
	472	Double_t AliEmcalJet::GetZ ( const AliVParticle* trk ) const
	473	{
	474	// Get Z of constituent trk
	475	return GetZ ( trk->Px(), trk->Py(), trk->Pz() );
	476	}
	477
f660c2d6	478	//__________________________________________________________________________________________________
	479	AliVParticle* AliEmcalJet::GetLeadingTrack(TClonesArray *tracks) const
	480	{
	481	AliVParticle* maxTrack = 0;
	482	for (Int_t i = 0; i < GetNumberOfTracks(); i++) {
	483	AliVParticle *track = TrackAt(i, tracks);
59ebe296	484	if (!track) {
3c9775d9	485	AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)",
3c9775d9	486	i,tracks->GetName(),TrackAt(i),tracks->GetEntriesFast()));
59ebe296	487	continue;
59ebe296	488	}
f660c2d6	489	if (!maxTrack \|\| track->Pt() > maxTrack->Pt())
	490	maxTrack = track;
	491	}
	492
	493	return maxTrack;
	494	}
	495
	496	//__________________________________________________________________________________________________
	497	AliVCluster* AliEmcalJet::GetLeadingCluster(TClonesArray *clusters) const
	498	{
	499	AliVCluster* maxCluster = 0;
	500	for (Int_t i = 0; i < GetNumberOfClusters(); i++) {
	501	AliVCluster *cluster = ClusterAt(i, clusters);
59ebe296	502	if (!cluster) {
3c9775d9	503	AliError(Form("Unable to find jet cluster %d in collection %s (pos in collection %d, max %d)",
3c9775d9	504	i,clusters->GetName(),ClusterAt(i),clusters->GetEntriesFast()));
59ebe296	505	continue;
59ebe296	506	}
f660c2d6	507	if (!maxCluster \|\| cluster->E() > maxCluster->E())
	508	maxCluster = cluster;
	509	}
	510
	511	return maxCluster;
	512	}
2103dc6a	513
	514	//__________________________________________________________________________________________________
	515	void AliEmcalJet::ResetMatching()
	516	{
	517	fClosestJets[0] = 0;
	518	fClosestJets[1] = 0;
	519	fClosestJetsDist[0] = 999;
	520	fClosestJetsDist[1] = 999;
	521	fMatched = 2;
	522	}
726c9488	523
	524	//__________________________________________________________________________________________________
	525	void AliEmcalJet::PrintGR() {
	526	for(Int_t i = 0; i<fGRNumerator.GetSize(); i++) {
	527	Printf("num[%d] = %f",i,fGRNumerator.At(i));
	528	}
	529	}