[u/mrichter/AliRoot.git] / PWGGA / 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;
  }

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

  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 (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 (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 (fClus) {
    if (fIsEmcPart) {
      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);  
      }
    } else { /*CaloClusters given as input*/
      const Int_t Nclus = fClus->GetEntries();
      for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
        AliVCluster *c = static_cast<AliVCluster*>(fClus->At(iClus));
        if (!c)
          continue;
        if (!c->IsEMCAL())
          continue;
        TLorentzVector nPart;
        c->GetMomentum(nPart, vertex);
        if (nPart.Pt() < fMinJetClusPt) 
          continue;
        // offset of 100 to skip ghost particles uid = -1
        fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);  
      }
    }
  }

  // run jet finder
  fjw.Run();

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

  // 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());
    jet->SetArea(fjw.GetJetArea(ij));
    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++;
  }
}

//________________________________________________________________________
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
  if (!(fEvent->FindListObject(fJetsName)))
    fEvent->AddObject(fJets);
  else {
    AliError(Form("%s: Object with name %s already in event! Returning", fJetsName.Data(), GetName()));
    return 0;
  }

  if ((fType == 0) || (fType == 1)) {
    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;
      }
    }
    TString objname(fTracks->GetClass()->GetName());
    TClass cls(objname);
    if (cls.InheritsFrom("AliMCParticle"))
      fIsMcPart = 1;
  }

  if ((fType == 0) || (fType == 2)) {
    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;
      }
    }
    TString objname(fClus->GetClass()->GetName());
    TClass cls(objname);
    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
48d874ff	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	}
	56
	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
63206072	113	// delete jet output
	114	fJets->Delete();
	115
04905231	116	FindJets();
48d874ff	117	}
	118
	119	//________________________________________________________________________
914d486c	120	void AliEmcalJetTask::Terminate(Option_t *)
48d874ff	121	{
48d874ff	122	// Called once at the end of the analysis.
48d874ff	123	}
	124
	125	//________________________________________________________________________
04905231	126	void AliEmcalJetTask::FindJets()
48d874ff	127	{
	128	// Find jets.
	129
04905231	130	if (!fTracks && !fClus)
d03084cd	131	return;
d03084cd	132
48d874ff	133	TString name("kt");
48d874ff	134	fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
04905231	135	if (fAlgo>=1) {
48d874ff	136	name = "antikt";
	137	jalgo = fastjet::antikt_algorithm;
	138	}
	139
04905231	140	// setup fj wrapper
48d874ff	141	AliFJWrapper fjw(name, name);
f5f3c8e9	142	fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
04905231	143	fjw.SetR(fRadius);
6ea168d0	144	fjw.SetAlgorithm(jalgo);
04905231	145	fjw.SetMaxRap(1);
48d874ff	146	fjw.Clear();
48d874ff	147
04905231	148	// get primary vertex
	149	Double_t vertex[3] = {0, 0, 0};
	150	fEvent->GetPrimaryVertex()->GetXYZ(vertex);
	151
	152	if (fTracks) {
	153	const Int_t Ntracks = fTracks->GetEntries();
48d874ff	154	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
04905231	155	AliVParticle t = static_cast<AliVParticle>(fTracks->At(iTracks));
48d874ff	156	if (!t)
48d874ff	157	continue;
1772fe65	158	if (t->Pt() < fMinJetTrackPt)
914d486c	159	continue;
04905231	160	Double_t eta = t->Eta();
	161	Double_t phi = t->Phi();
	162	if ((eta<fEtaMin) && (eta>fEtaMax) &&
	163	(phi<fPhiMin) && (phi<fPhiMax))
	164	continue;
a3c8c8c8	165	// offset of 100 for consistency with cluster ids
a3c8c8c8	166	fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
48d874ff	167	}
48d874ff	168	}
6ea168d0	169
04905231	170	if (fClus) {
	171	if (fIsEmcPart) {
	172	const Int_t Nclus = fClus->GetEntries();
	173	for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
	174	AliVCluster *c = 0;
	175	Double_t cEta=0,cPhi=0,cPt=0;
	176	Double_t cPx=0,cPy=0,cPz=0;
	177	if (fIsEmcPart) {
	178	AliEmcalParticle ep = static_cast<AliEmcalParticle>(fClus->At(iClus));
	179	if (!ep)
	180	continue;
	181	c = ep->GetCluster();
	182	if (!c)
	183	continue;
	184	cEta = ep->Eta();
	185	cPhi = ep->Phi();
	186	cPt = ep->Pt();
	187	cPx = ep->Px();
	188	cPy = ep->Py();
	189	cPz = ep->Pz();
	190	} else {
	191	c = static_cast<AliVCluster*>(fClus->At(iClus));
	192	if (!c)
	193	continue;
	194	TLorentzVector nP;
	195	c->GetMomentum(nP, vertex);
	196	cEta = nP.Eta();
	197	cPhi = nP.Phi();
	198	cPt = nP.Pt();
	199	cPx = nP.Px();
	200	cPy = nP.Py();
	201	cPz = nP.Pz();
	202	}
	203	if (!c->IsEMCAL())
	204	continue;
	205	if (cPt < fMinJetClusPt)
	206	continue;
	207	if ((cEta<fEtaMin) && (cEta>fEtaMax) &&
	208	(cPhi<fPhiMin) && (cPhi<fPhiMax))
	209	continue;
	210	// offset of 100 to skip ghost particles uid = -1
	211	fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPxcPx+cPycPy+cPz*cPz), -iClus - 100);
	212	}
	213	} else { /CaloClusters given as input/
	214	const Int_t Nclus = fClus->GetEntries();
	215	for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
	216	AliVCluster c = static_cast<AliVCluster>(fClus->At(iClus));
	217	if (!c)
	218	continue;
	219	if (!c->IsEMCAL())
	220	continue;
	221	TLorentzVector nPart;
	222	c->GetMomentum(nPart, vertex);
	223	if (nPart.Pt() < fMinJetClusPt)
	224	continue;
	225	// offset of 100 to skip ghost particles uid = -1
	226	fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);
	227	}
48d874ff	228	}
48d874ff	229	}
04905231	230
7efbea04	231	// run jet finder
48d874ff	232	fjw.Run();
f5f3c8e9	233
	234	// get geometry
	235	AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
	236	if (!geom) {
02c7e943	237	AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9	238	return;
f5f3c8e9	239	}
04905231	240
04905231	241	// loop over fastjet jets
48d874ff	242	std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
b65fac7a	243	for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9	244	if (jets_incl[ij].perp()<fMinJetPt)
	245	continue;
	246	if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff	247	continue;
04905231	248
914d486c	249	AliEmcalJet jet = new ((fJets)[jetCount])
914d486c	250	AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
04905231	251
04905231	252	// loop over constituents
c64cb1f6	253	std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
04905231	254	jet->SetNumberOfTracks(constituents.size());
	255	jet->SetNumberOfClusters(constituents.size());
	256
b65fac7a	257	Int_t nt = 0;
b65fac7a	258	Int_t nc = 0;
f5f3c8e9	259	Double_t neutralE = 0;
111318e8	260	Double_t maxCh = 0;
111318e8	261	Double_t maxNe = 0;
f5f3c8e9	262	Int_t gall = 0;
f5f3c8e9	263	Int_t gemc = 0;
04905231	264	Int_t cemc = 0;
111318e8	265	Int_t ncharged = 0;
111318e8	266	Int_t nneutral = 0;
02c7e943	267	Double_t mcpt = 0;
04905231	268	Double_t emcpt = 0;
b65fac7a	269
b65fac7a	270	for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff	271	Int_t uid = constituents[ic].user_index();
b65fac7a	272
04905231	273	if ((uid == -1) /&& (constituents[ic].kt2() < 1e-25)/) { //ghost particle
f5f3c8e9	274	++gall;
04905231	275	Double_t gphi = constituents[ic].phi();
	276	if (gphi<0)
	277	gphi += TMath::TwoPi();
	278	gphi *= TMath::RadToDeg();
f5f3c8e9	279	Double_t geta = constituents[ic].eta();
b65fac7a	280	if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
b65fac7a	281	(geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
04905231	282	++gemc;
04905231	283	} else if ((uid > 0) && fTracks) { // track constituent
111318e8	284	Int_t tid = uid - 100;
04905231	285	AliVParticle t = static_cast<AliVParticle>(fTracks->At(tid));
	286	if (!t)
	287	continue;
	288	Double_t cEta = t->Eta();
	289	Double_t cPhi = t->Phi();
	290	Double_t cPt = t->Pt();
	291	Double_t cP = t->P();
	292	if (t->Charge() == 0) {
	293	neutralE += cP;
	294	++nneutral;
	295	if (cPt > maxNe)
	296	maxNe = cPt;
	297	} else {
	298	++ncharged;
	299	if (cPt > maxCh)
	300	maxCh = cPt;
	301	}
	302
	303	if (fIsMcPart \|\| t->GetLabel() == 100) // check if MC particle
	304	mcpt += cPt;
	305
	306	if (cPhi<0)
	307	cPhi += TMath::TwoPi();
	308	cPhi *= TMath::RadToDeg();
	309	if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
	310	(cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
	311	emcpt += cPt;
	312	++cemc;
d03084cd	313	}
04905231	314
	315	jet->AddTrackAt(tid, nt);
	316	++nt;
	317	} else if (fClus) { // cluster constituent
111318e8	318	Int_t cid = -uid - 100;
04905231	319	AliVCluster *c = 0;
	320	Double_t cEta=0,cPhi=0,cPt=0,cP=0;
	321	if (fIsEmcPart) {
	322	AliEmcalParticle ep = static_cast<AliEmcalParticle>(fClus->At(cid));
	323	if (!ep)
	324	continue;
	325	c = ep->GetCluster();
	326	if (!c)
	327	continue;
	328	cEta = ep->Eta();
	329	cPhi = ep->Phi();
	330	cPt = ep->Pt();
	331	cP = ep->P();
	332	} else {
	333	c = static_cast<AliVCluster*>(fClus->At(cid));
	334	if (!c)
	335	continue;
	336	TLorentzVector nP;
	337	c->GetMomentum(nP, vertex);
	338	cEta = nP.Eta();
	339	cPhi = nP.Phi();
	340	cPt = nP.Pt();
	341	cP = nP.P();
	342	}
	343
	344	neutralE += cP;
	345	if (cPt > maxNe)
	346	maxNe = cPt;
	347
	348	if (c->Chi2() == 100) // MC particle
	349	mcpt += cPt;
	350
	351	if (cPhi<0)
	352	cPhi += TMath::TwoPi();
	353	cPhi *= TMath::RadToDeg();
	354	if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
	355	(cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
	356	emcpt += cPt;
	357	++cemc;
	358	}
	359
	360	jet->AddClusterAt(cid, nc);
	361	++nc;
	362	++nneutral;
02c7e943	363	} else {
	364	AliError(Form("%s: No logical way to end up here.", GetName()));
	365	continue;
48d874ff	366	}
04905231	367	} /* end of constituent loop */
04905231	368
35789a2d	369	jet->SetNumberOfTracks(nt);
35789a2d	370	jet->SetNumberOfClusters(nc);
f5f3c8e9	371	jet->SortConstituents();
111318e8	372	jet->SetMaxNeutralPt(maxNe);
111318e8	373	jet->SetMaxChargedPt(maxCh);
b65fac7a	374	jet->SetNEF(neutralE / jet->E());
f5f3c8e9	375	jet->SetArea(fjw.GetJetArea(ij));
111318e8	376	jet->SetNumberOfCharged(ncharged);
111318e8	377	jet->SetNumberOfNeutrals(nneutral);
02c7e943	378	jet->SetMCPt(mcpt);
04905231	379	jet->SetNEmc(cemc);
	380	jet->SetPtEmc(emcpt);
	381
b65fac7a	382	if (gall > 0)
b65fac7a	383	jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9	384	else
f5f3c8e9	385	jet->SetAreaEmc(-1);
b65fac7a	386	if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65	387	(jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
	388	(jet->Eta() > geom->GetArm1EtaMin()) &&
	389	(jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9	390	jet->SetAxisInEmcal(kTRUE);
48d874ff	391	jetCount++;
	392	}
	393	}
04905231	394
	395	//________________________________________________________________________
	396	Bool_t AliEmcalJetTask::DoInit()
	397	{
	398	// Init. Return true if successful.
	399
	400	// get input collections
	401	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	402
	403	// get the event
	404	fEvent = InputEvent();
	405	if (!fEvent) {
	406	AliError(Form("%s: Could not retrieve event! Returning", GetName()));
	407	return 0;
	408	}
	409
	410	// add jets to event if not yet there
	411	if (!(fEvent->FindListObject(fJetsName)))
	412	fEvent->AddObject(fJets);
	413	else {
	414	AliError(Form("%s: Object with name %s already in event! Returning", fJetsName.Data(), GetName()));
	415	return 0;
	416	}
	417
	418	if ((fType == 0) \|\| (fType == 1)) {
	419	if (fTracksName == "Tracks")
	420	am->LoadBranch("Tracks");
	421	if (!fTracks && !fTracksName.IsNull()) {
	422	fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
	423	if (!fTracks) {
	424	AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
	425	return 0;
	426	}
	427	}
	428	TString objname(fTracks->GetClass()->GetName());
	429	TClass cls(objname);
	430	if (cls.InheritsFrom("AliMCParticle"))
	431	fIsMcPart = 1;
	432	}
	433
	434	if ((fType == 0) \|\| (fType == 2)) {
	435	if (fCaloName == "CaloClusters")
	436	am->LoadBranch("CaloClusters");
	437	if (!fClus && !fCaloName.IsNull()) {
	438	fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
	439	if (!fClus) {
	440	AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
	441	return 0;
	442	}
	443	}
	444	TString objname(fClus->GetClass()->GetName());
	445	TClass cls(objname);
	446	if (cls.InheritsFrom("AliEmcalParticle"))
	447	fIsEmcPart = 1;
	448	}
	449
	450	return 1;
	451	}