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

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

#include <vector>
#include "AliEmcalJetTask.h"

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

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

ClassImp(AliEmcalJetTask)

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask() : 
  AliAnalysisTaskSE("AliEmcalJetTask"),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fAlgo(1),
  fRadius(0.4),
  fType(0),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fPhiMin(0),
  fPhiMax(TMath::TwoPi()),
  fEtaMin(-1),
  fEtaMax(1),
  fMinJetArea(0.01),
  fMinJetPt(1.0),
  fIsInit(0),
  fIsMcPart(0),
  fIsEmcPart(0),
  fJets(0),
  fEvent(0),
  fTracks(0),
  fClus(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask(const char *name) : 
  AliAnalysisTaskSE(name),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fAlgo(1),
  fRadius(0.4),
  fType(0),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fPhiMin(0),
  fPhiMax(TMath::TwoPi()),
  fEtaMin(-1),
  fEtaMax(1),
  fMinJetArea(0.01),
  fMinJetPt(1.0),
  fIsInit(0),
  fIsMcPart(0),
  fIsEmcPart(0),
  fJets(0),
  fEvent(0),
  fTracks(0),
  fClus(0)
{
  // Standard constructor.

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

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

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

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

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

  if (!fIsInit) {
    if (!DoInit())
      return;
    fIsInit = kTRUE;
  }

  FindJets();
}

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

//________________________________________________________________________
void AliEmcalJetTask::FindJets()
{
  // Find jets.

  if (!fTracks && !fClus)
    return;

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

  // setup fj wrapper
  AliFJWrapper fjw(name, name);
  fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
  fjw.SetR(fRadius);
  fjw.SetAlgorithm(jalgo);  
  fjw.SetMaxRap(1);
  fjw.Clear();

  // get primary vertex
  Double_t vertex[3] = {0, 0, 0};
  fEvent->GetPrimaryVertex()->GetXYZ(vertex);

  if ((fIsMcPart || fType == 0 || fType == 1) && fTracks) {
    const Int_t Ntracks = fTracks->GetEntries();
    for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
      AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
      if (!t)
        continue;
      if (fIsMcPart && fType == 1 && t->Charge() == 0)
	continue;
      if (fIsMcPart && fType == 2 && t->Charge() != 0)
	continue;
      if (t->Pt() < fMinJetTrackPt) 
        continue;
      Double_t eta = t->Eta();
      Double_t phi = t->Phi();
      if ((eta<fEtaMin) && (eta>fEtaMax) &&
          (phi<fPhiMin) && (phi<fPhiMax))
        continue;

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

  if ((fType == 0 || fType == 2) && fClus) {
    const Int_t Nclus = fClus->GetEntries();
    for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
      AliVCluster *c = 0;
      Double_t cEta=0,cPhi=0,cPt=0;
      Double_t cPx=0,cPy=0,cPz=0;
      if (fIsEmcPart) {
	AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
	if (!ep)
	  continue;
	c = ep->GetCluster();
	if (!c)
            continue;
	cEta = ep->Eta();
	cPhi = ep->Phi();
	cPt  = ep->Pt();
	cPx  = ep->Px();
	cPy  = ep->Py();
	cPz  = ep->Pz();
      } else {
	c = static_cast<AliVCluster*>(fClus->At(iClus));
	if (!c)
	  continue;
	TLorentzVector nP;
	c->GetMomentum(nP, vertex);
	cEta = nP.Eta();
	cPhi = nP.Phi();
	cPt  = nP.Pt();
	cPx  = nP.Px();
	cPy  = nP.Py();
	cPz  = nP.Pz();
      }
      if (!c->IsEMCAL())
	continue;
      if (cPt < fMinJetClusPt) 
	continue;
      if ((cEta<fEtaMin) && (cEta>fEtaMax) &&
	  (cPhi<fPhiMin) && (cPhi<fPhiMax))
	continue;
      // offset of 100 to skip ghost particles uid = -1
      fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);   
    }
  }
  
  // run jet finder
  fjw.Run();

  // get geometry
  AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
  if (!geom) {
    AliFatal(Form("%s: Can not create geometry", GetName()));
    return;
  }

  // loop over fastjet jets
  std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
  for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
    if (jets_incl[ij].perp()<fMinJetPt) 
      continue;
    if (fjw.GetJetArea(ij)<fMinJetArea)
      continue;

    AliEmcalJet *jet = new ((*fJets)[jetCount]) 
      AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());

    // loop over constituents
    std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
    jet->SetNumberOfTracks(constituents.size());
    jet->SetNumberOfClusters(constituents.size());

    Int_t nt            = 0;
    Int_t nc            = 0;
    Double_t neutralE   = 0;
    Double_t maxCh      = 0;
    Double_t maxNe      = 0;
    Int_t gall          = 0;
    Int_t gemc          = 0;
    Int_t cemc          = 0;
    Int_t ncharged      = 0;
    Int_t nneutral      = 0;
    Double_t mcpt       = 0;
    Double_t emcpt      = 0;

    for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
      Int_t uid = constituents[ic].user_index();

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

        if (fIsMcPart || t->GetLabel() == 100) // check if MC particle
          mcpt += cPt;

        if (cPhi<0) 
          cPhi += TMath::TwoPi();
        cPhi *= TMath::RadToDeg();
        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
          emcpt += cPt;
          ++cemc;
        }

        jet->AddTrackAt(tid, nt);
        ++nt;
      } else if (fClus) { // cluster constituent
	Int_t cid = -uid - 100;
	AliVCluster *c = 0;
        Double_t cEta=0,cPhi=0,cPt=0,cP=0;
        if (fIsEmcPart) {
          AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
          if (!ep)
            continue;
          c = ep->GetCluster();
          if (!c)
            continue;
          cEta = ep->Eta();
          cPhi = ep->Phi();
          cPt  = ep->Pt();
          cP   = ep->P();
        } else {
          c = static_cast<AliVCluster*>(fClus->At(cid));
          if (!c)
            continue;
          TLorentzVector nP;
          c->GetMomentum(nP, vertex);
          cEta = nP.Eta();
          cPhi = nP.Phi();
          cPt  = nP.Pt();
          cP   = nP.P();
        }

        neutralE += cP;
        if (cPt > maxNe)
          maxNe = cPt;

        if (c->Chi2() == 100) // MC particle
          mcpt += cPt;

        if (cPhi<0) 
          cPhi += TMath::TwoPi();
        cPhi *= TMath::RadToDeg();
        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
          emcpt += cPt;
          ++cemc;
        }

        jet->AddClusterAt(cid, nc);
        ++nc;
        ++nneutral;
      } else {
        AliError(Form("%s: No logical way to end up here.", GetName()));
        continue;
      }
    } /* end of constituent loop */

    jet->SetNumberOfTracks(nt);
    jet->SetNumberOfClusters(nc);
    jet->SortConstituents();
    jet->SetMaxNeutralPt(maxNe);
    jet->SetMaxChargedPt(maxCh);
    jet->SetNEF(neutralE / jet->E());
    fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
    jet->SetArea(area.perp());
    jet->SetAreaEta(area.eta());
    jet->SetAreaPhi(area.phi());
    jet->SetNumberOfCharged(ncharged);
    jet->SetNumberOfNeutrals(nneutral);
    jet->SetMCPt(mcpt);
    jet->SetNEmc(cemc);
    jet->SetPtEmc(emcpt);

    if (gall > 0)
      jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
    else 
      jet->SetAreaEmc(-1);
    if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && 
	(jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
	(jet->Eta() > geom->GetArm1EtaMin()) && 
	(jet->Eta() < geom->GetArm1EtaMax()))
      jet->SetAxisInEmcal(kTRUE);
    jetCount++;
  }
  fJets->Sort();
}

//________________________________________________________________________
Bool_t AliEmcalJetTask::DoInit()
{
  // Init. Return true if successful.

  // get input collections
  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();

  // get the event
  fEvent = InputEvent();
  if (!fEvent) {
    AliError(Form("%s: Could not retrieve event! Returning", GetName()));
    return 0;
  }

  // add jets to event if not yet there
  fJets->Delete();
  if (!(fEvent->FindListObject(fJetsName)))
    fEvent->AddObject(fJets);
  else {
    AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
    return 0;
  }

  if (fTracksName == "Tracks")
    am->LoadBranch("Tracks");
  if (!fTracks && !fTracksName.IsNull()) {
    fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
    if (!fTracks) {
      AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
      return 0;
    }
  }
  if (fTracks) {
    TClass cls(fTracks->GetClass()->GetName());
    if (cls.InheritsFrom("AliMCParticle"))
      fIsMcPart = 1;
  }
  
  if (fCaloName == "CaloClusters")
    am->LoadBranch("CaloClusters");
  if (!fClus && !fCaloName.IsNull()) {
    fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
    if (!fClus) {
      AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
      return 0;
    }
  }
  if (fClus) {
    TClass cls(fClus->GetClass()->GetName());
    if (cls.InheritsFrom("AliEmcalParticle"))
      fIsEmcPart = 1;
  }

  return 1;
}
Commit	Line	Data
48d874ff	1	// $Id$
542ac102	2	//
	3	// Emcal jet finder task.
	4	//
b65fac7a	5	// Authors: C.Loizides, S.Aiola
48d874ff	6
c64cb1f6	7	#include <vector>
04905231	8	#include "AliEmcalJetTask.h"
04905231	9
6ea168d0	10	#include <TChain.h>
48d874ff	11	#include <TClonesArray.h>
6ea168d0	12	#include <TList.h>
6ea168d0	13	#include <TLorentzVector.h>
b65fac7a	14
48d874ff	15	#include "AliAnalysisManager.h"
6ea168d0	16	#include "AliCentrality.h"
f5f3c8e9	17	#include "AliEMCALGeometry.h"
04905231	18	#include "AliESDEvent.h"
f5f3c8e9	19	#include "AliEmcalJet.h"
04905231	20	#include "AliEmcalParticle.h"
63206072	21	#include "AliFJWrapper.h"
04905231	22	#include "AliMCEvent.h"
35789a2d	23	#include "AliVCluster.h"
b65fac7a	24	#include "AliVEvent.h"
04905231	25	#include "AliVParticle.h"
48d874ff	26
914d486c	27	ClassImp(AliEmcalJetTask)
48d874ff	28
6ea168d0	29	//________________________________________________________________________
914d486c	30	AliEmcalJetTask::AliEmcalJetTask() :
914d486c	31	AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0	32	fTracksName("Tracks"),
	33	fCaloName("CaloClusters"),
	34	fJetsName("Jets"),
	35	fAlgo(1),
	36	fRadius(0.4),
	37	fType(0),
914d486c	38	fMinJetTrackPt(0.15),
914d486c	39	fMinJetClusPt(0.15),
04905231	40	fPhiMin(0),
	41	fPhiMax(TMath::TwoPi()),
	42	fEtaMin(-1),
	43	fEtaMax(1),
f5f3c8e9	44	fMinJetArea(0.01),
f5f3c8e9	45	fMinJetPt(1.0),
04905231	46	fIsInit(0),
	47	fIsMcPart(0),
	48	fIsEmcPart(0),
b65fac7a	49	fJets(0),
04905231	50	fEvent(0),
	51	fTracks(0),
	52	fClus(0)
6ea168d0	53	{
6ea168d0	54	// Default constructor.
6ea168d0	55	}
6ea168d0	56
48d874ff	57	//________________________________________________________________________
914d486c	58	AliEmcalJetTask::AliEmcalJetTask(const char *name) :
63206072	59	AliAnalysisTaskSE(name),
7efbea04	60	fTracksName("Tracks"),
48d874ff	61	fCaloName("CaloClusters"),
7efbea04	62	fJetsName("Jets"),
48d874ff	63	fAlgo(1),
	64	fRadius(0.4),
	65	fType(0),
6ea168d0	66	fMinJetTrackPt(0.15),
914d486c	67	fMinJetClusPt(0.15),
04905231	68	fPhiMin(0),
	69	fPhiMax(TMath::TwoPi()),
	70	fEtaMin(-1),
	71	fEtaMax(1),
f5f3c8e9	72	fMinJetArea(0.01),
f5f3c8e9	73	fMinJetPt(1.0),
04905231	74	fIsInit(0),
	75	fIsMcPart(0),
	76	fIsEmcPart(0),
b65fac7a	77	fJets(0),
04905231	78	fEvent(0),
	79	fTracks(0),
	80	fClus(0)
48d874ff	81	{
48d874ff	82	// Standard constructor.
7efbea04	83
7efbea04	84	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff	85	}
	86
	87	//________________________________________________________________________
914d486c	88	AliEmcalJetTask::~AliEmcalJetTask()
48d874ff	89	{
	90	// Destructor
	91	}
	92
	93	//________________________________________________________________________
914d486c	94	void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff	95	{
	96	// Create user objects.
	97
914d486c	98	fJets = new TClonesArray("AliEmcalJet");
48d874ff	99	fJets->SetName(fJetsName);
48d874ff	100	}
	101
	102	//________________________________________________________________________
914d486c	103	void AliEmcalJetTask::UserExec(Option_t *)
48d874ff	104	{
	105	// Main loop, called for each event.
	106
04905231	107	if (!fIsInit) {
	108	if (!DoInit())
	109	return;
	110	fIsInit = kTRUE;
1772fe65	111	}
b65fac7a	112
04905231	113	FindJets();
48d874ff	114	}
	115
	116	//________________________________________________________________________
914d486c	117	void AliEmcalJetTask::Terminate(Option_t *)
48d874ff	118	{
48d874ff	119	// Called once at the end of the analysis.
48d874ff	120	}
	121
	122	//________________________________________________________________________
04905231	123	void AliEmcalJetTask::FindJets()
48d874ff	124	{
	125	// Find jets.
	126
04905231	127	if (!fTracks && !fClus)
d03084cd	128	return;
d03084cd	129
48d874ff	130	TString name("kt");
48d874ff	131	fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
04905231	132	if (fAlgo>=1) {
48d874ff	133	name = "antikt";
	134	jalgo = fastjet::antikt_algorithm;
	135	}
	136
04905231	137	// setup fj wrapper
48d874ff	138	AliFJWrapper fjw(name, name);
f5f3c8e9	139	fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
04905231	140	fjw.SetR(fRadius);
6ea168d0	141	fjw.SetAlgorithm(jalgo);
04905231	142	fjw.SetMaxRap(1);
48d874ff	143	fjw.Clear();
48d874ff	144
04905231	145	// get primary vertex
	146	Double_t vertex[3] = {0, 0, 0};
	147	fEvent->GetPrimaryVertex()->GetXYZ(vertex);
	148
da67b88a	149	if ((fIsMcPart \|\| fType == 0 \|\| fType == 1) && fTracks) {
04905231	150	const Int_t Ntracks = fTracks->GetEntries();
48d874ff	151	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
04905231	152	AliVParticle t = static_cast<AliVParticle>(fTracks->At(iTracks));
48d874ff	153	if (!t)
48d874ff	154	continue;
5d950148	155	if (fIsMcPart && fType == 1 && t->Charge() == 0)
da67b88a	156	continue;
5d950148	157	if (fIsMcPart && fType == 2 && t->Charge() != 0)
da67b88a	158	continue;
1772fe65	159	if (t->Pt() < fMinJetTrackPt)
914d486c	160	continue;
04905231	161	Double_t eta = t->Eta();
	162	Double_t phi = t->Phi();
	163	if ((eta<fEtaMin) && (eta>fEtaMax) &&
	164	(phi<fPhiMin) && (phi<fPhiMax))
	165	continue;
da67b88a	166
a3c8c8c8	167	// offset of 100 for consistency with cluster ids
a3c8c8c8	168	fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
48d874ff	169	}
48d874ff	170	}
6ea168d0	171
da67b88a	172	if ((fType == 0 \|\| fType == 2) && fClus) {
	173	const Int_t Nclus = fClus->GetEntries();
	174	for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
	175	AliVCluster *c = 0;
	176	Double_t cEta=0,cPhi=0,cPt=0;
	177	Double_t cPx=0,cPy=0,cPz=0;
	178	if (fIsEmcPart) {
	179	AliEmcalParticle ep = static_cast<AliEmcalParticle>(fClus->At(iClus));
	180	if (!ep)
	181	continue;
	182	c = ep->GetCluster();
	183	if (!c)
04905231	184	continue;
da67b88a	185	cEta = ep->Eta();
	186	cPhi = ep->Phi();
	187	cPt = ep->Pt();
	188	cPx = ep->Px();
	189	cPy = ep->Py();
	190	cPz = ep->Pz();
	191	} else {
	192	c = static_cast<AliVCluster*>(fClus->At(iClus));
	193	if (!c)
	194	continue;
	195	TLorentzVector nP;
	196	c->GetMomentum(nP, vertex);
	197	cEta = nP.Eta();
	198	cPhi = nP.Phi();
	199	cPt = nP.Pt();
	200	cPx = nP.Px();
	201	cPy = nP.Py();
	202	cPz = nP.Pz();
04905231	203	}
da67b88a	204	if (!c->IsEMCAL())
	205	continue;
	206	if (cPt < fMinJetClusPt)
	207	continue;
	208	if ((cEta<fEtaMin) && (cEta>fEtaMax) &&
	209	(cPhi<fPhiMin) && (cPhi<fPhiMax))
	210	continue;
	211	// offset of 100 to skip ghost particles uid = -1
	212	fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPxcPx+cPycPy+cPz*cPz), -iClus - 100);
48d874ff	213	}
48d874ff	214	}
da67b88a	215
7efbea04	216	// run jet finder
48d874ff	217	fjw.Run();
f5f3c8e9	218
	219	// get geometry
	220	AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
	221	if (!geom) {
02c7e943	222	AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9	223	return;
f5f3c8e9	224	}
04905231	225
04905231	226	// loop over fastjet jets
48d874ff	227	std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
b65fac7a	228	for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9	229	if (jets_incl[ij].perp()<fMinJetPt)
	230	continue;
	231	if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff	232	continue;
04905231	233
914d486c	234	AliEmcalJet jet = new ((fJets)[jetCount])
914d486c	235	AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
04905231	236
04905231	237	// loop over constituents
c64cb1f6	238	std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
04905231	239	jet->SetNumberOfTracks(constituents.size());
	240	jet->SetNumberOfClusters(constituents.size());
	241
b65fac7a	242	Int_t nt = 0;
b65fac7a	243	Int_t nc = 0;
f5f3c8e9	244	Double_t neutralE = 0;
111318e8	245	Double_t maxCh = 0;
111318e8	246	Double_t maxNe = 0;
f5f3c8e9	247	Int_t gall = 0;
f5f3c8e9	248	Int_t gemc = 0;
04905231	249	Int_t cemc = 0;
111318e8	250	Int_t ncharged = 0;
111318e8	251	Int_t nneutral = 0;
02c7e943	252	Double_t mcpt = 0;
04905231	253	Double_t emcpt = 0;
b65fac7a	254
b65fac7a	255	for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff	256	Int_t uid = constituents[ic].user_index();
b65fac7a	257
04905231	258	if ((uid == -1) /&& (constituents[ic].kt2() < 1e-25)/) { //ghost particle
f5f3c8e9	259	++gall;
04905231	260	Double_t gphi = constituents[ic].phi();
	261	if (gphi<0)
	262	gphi += TMath::TwoPi();
	263	gphi *= TMath::RadToDeg();
f5f3c8e9	264	Double_t geta = constituents[ic].eta();
b65fac7a	265	if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
b65fac7a	266	(geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
04905231	267	++gemc;
04905231	268	} else if ((uid > 0) && fTracks) { // track constituent
111318e8	269	Int_t tid = uid - 100;
04905231	270	AliVParticle t = static_cast<AliVParticle>(fTracks->At(tid));
	271	if (!t)
	272	continue;
	273	Double_t cEta = t->Eta();
	274	Double_t cPhi = t->Phi();
	275	Double_t cPt = t->Pt();
	276	Double_t cP = t->P();
	277	if (t->Charge() == 0) {
	278	neutralE += cP;
	279	++nneutral;
	280	if (cPt > maxNe)
	281	maxNe = cPt;
	282	} else {
	283	++ncharged;
	284	if (cPt > maxCh)
	285	maxCh = cPt;
	286	}
	287
	288	if (fIsMcPart \|\| t->GetLabel() == 100) // check if MC particle
	289	mcpt += cPt;
	290
	291	if (cPhi<0)
	292	cPhi += TMath::TwoPi();
	293	cPhi *= TMath::RadToDeg();
	294	if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
	295	(cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
	296	emcpt += cPt;
	297	++cemc;
d03084cd	298	}
04905231	299
	300	jet->AddTrackAt(tid, nt);
	301	++nt;
	302	} else if (fClus) { // cluster constituent
111318e8	303	Int_t cid = -uid - 100;
04905231	304	AliVCluster *c = 0;
	305	Double_t cEta=0,cPhi=0,cPt=0,cP=0;
	306	if (fIsEmcPart) {
	307	AliEmcalParticle ep = static_cast<AliEmcalParticle>(fClus->At(cid));
	308	if (!ep)
	309	continue;
	310	c = ep->GetCluster();
	311	if (!c)
	312	continue;
	313	cEta = ep->Eta();
	314	cPhi = ep->Phi();
	315	cPt = ep->Pt();
	316	cP = ep->P();
	317	} else {
	318	c = static_cast<AliVCluster*>(fClus->At(cid));
	319	if (!c)
	320	continue;
	321	TLorentzVector nP;
	322	c->GetMomentum(nP, vertex);
	323	cEta = nP.Eta();
	324	cPhi = nP.Phi();
	325	cPt = nP.Pt();
	326	cP = nP.P();
	327	}
	328
	329	neutralE += cP;
	330	if (cPt > maxNe)
	331	maxNe = cPt;
	332
	333	if (c->Chi2() == 100) // MC particle
	334	mcpt += cPt;
	335
	336	if (cPhi<0)
	337	cPhi += TMath::TwoPi();
	338	cPhi *= TMath::RadToDeg();
	339	if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
	340	(cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
	341	emcpt += cPt;
	342	++cemc;
	343	}
	344
	345	jet->AddClusterAt(cid, nc);
	346	++nc;
	347	++nneutral;
02c7e943	348	} else {
	349	AliError(Form("%s: No logical way to end up here.", GetName()));
	350	continue;
48d874ff	351	}
04905231	352	} /* end of constituent loop */
04905231	353
35789a2d	354	jet->SetNumberOfTracks(nt);
35789a2d	355	jet->SetNumberOfClusters(nc);
f5f3c8e9	356	jet->SortConstituents();
111318e8	357	jet->SetMaxNeutralPt(maxNe);
111318e8	358	jet->SetMaxChargedPt(maxCh);
b65fac7a	359	jet->SetNEF(neutralE / jet->E());
a88b5c99	360	fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
	361	jet->SetArea(area.perp());
	362	jet->SetAreaEta(area.eta());
	363	jet->SetAreaPhi(area.phi());
111318e8	364	jet->SetNumberOfCharged(ncharged);
111318e8	365	jet->SetNumberOfNeutrals(nneutral);
02c7e943	366	jet->SetMCPt(mcpt);
04905231	367	jet->SetNEmc(cemc);
	368	jet->SetPtEmc(emcpt);
	369
b65fac7a	370	if (gall > 0)
b65fac7a	371	jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9	372	else
f5f3c8e9	373	jet->SetAreaEmc(-1);
b65fac7a	374	if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65	375	(jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
	376	(jet->Eta() > geom->GetArm1EtaMin()) &&
	377	(jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9	378	jet->SetAxisInEmcal(kTRUE);
48d874ff	379	jetCount++;
48d874ff	380	}
a88b5c99	381	fJets->Sort();
48d874ff	382	}
04905231	383
	384	//________________________________________________________________________
	385	Bool_t AliEmcalJetTask::DoInit()
	386	{
	387	// Init. Return true if successful.
	388
	389	// get input collections
	390	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	391
	392	// get the event
	393	fEvent = InputEvent();
	394	if (!fEvent) {
	395	AliError(Form("%s: Could not retrieve event! Returning", GetName()));
	396	return 0;
	397	}
	398
	399	// add jets to event if not yet there
71c1dd76	400	fJets->Delete();
04905231	401	if (!(fEvent->FindListObject(fJetsName)))
	402	fEvent->AddObject(fJets);
	403	else {
f4b49fa6	404	AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
04905231	405	return 0;
	406	}
	407
da67b88a	408	if (fTracksName == "Tracks")
	409	am->LoadBranch("Tracks");
	410	if (!fTracks && !fTracksName.IsNull()) {
	411	fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
	412	if (!fTracks) {
	413	AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
	414	return 0;
04905231	415	}
da67b88a	416	}
	417	if (fTracks) {
	418	TClass cls(fTracks->GetClass()->GetName());
04905231	419	if (cls.InheritsFrom("AliMCParticle"))
	420	fIsMcPart = 1;
	421	}
da67b88a	422
	423	if (fCaloName == "CaloClusters")
	424	am->LoadBranch("CaloClusters");
	425	if (!fClus && !fCaloName.IsNull()) {
	426	fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
	427	if (!fClus) {
	428	AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
	429	return 0;
04905231	430	}
da67b88a	431	}
	432	if (fClus) {
	433	TClass cls(fClus->GetClass()->GetName());
04905231	434	if (cls.InheritsFrom("AliEmcalParticle"))
	435	fIsEmcPart = 1;
	436	}
	437
	438	return 1;
	439	}