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

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

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

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

#include "AliEmcalJetTask.h"

ClassImp(AliEmcalJetTask)

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask() : 
  AliAnalysisTaskSE("AliEmcalJetTask"),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fMC(kFALSE),
  fAlgo(1),
  fRadius(0.4),
  fType(0),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fMinJetArea(0.01),
  fMinJetPt(1.0),
  fJets(0),
  fEvent(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliEmcalJetTask::AliEmcalJetTask(const char *name) : 
  AliAnalysisTaskSE(name),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fMC(kFALSE),
  fAlgo(1),
  fRadius(0.4),
  fType(0),
  fMinJetTrackPt(0.15),
  fMinJetClusPt(0.15),
  fMinJetArea(0.01),
  fMinJetPt(1.0),
  fJets(0),
  fEvent(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.

  // get the event
  if (fMC) 
    fEvent = MCEvent();
  else
    fEvent = InputEvent();

  // add jets to event if not yet there
  if (!(fEvent->FindListObject(fJetsName)))
    fEvent->AddObject(fJets);

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

  if (!fEvent) {
    AliError(Form("Could not get the event! fMC = %d", fMC));
    return;
  }

  // get input collections
  TClonesArray *tracks = 0;
  TClonesArray *clus   = 0;
  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();

  if ((fType == 0) || (fType == 1)) {
    if (fTracksName == "Tracks")
      am->LoadBranch("Tracks");
    tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
    if (!tracks) {
      AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
      return;
    }
  }
  if ((fType == 0) || (fType == 2)) {
    if (fCaloName == "CaloClusters")
      am->LoadBranch("CaloClusters");
    clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
    if (!clus) {
      AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
      return;
    }
  }

  // get centrality
  Float_t cent = -1; 
  AliCentrality *centrality = fEvent->GetCentrality();
  if (centrality)
    cent = centrality->GetCentralityPercentile("V0M");
  else
    cent=99; // probably pp data
  if (cent < 0) {
    AliError(Form("Centrality negative: %f", cent));
    return;
  }

  FindJets(tracks, clus, fAlgo, fRadius, cent);
}

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

//________________________________________________________________________
void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
{
  // Find jets.

  if (!tracks && !clus)
    return;

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

  AliFJWrapper fjw(name, name);
  fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
  fjw.SetR(radius);
  fjw.SetAlgorithm(jalgo);  
  fjw.SetMaxRap(0.9);
  fjw.Clear();

  if (tracks) {
    const Int_t Ntracks = tracks->GetEntries();
    for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
      AliVParticle *t = static_cast<AliVParticle*>(tracks->At(iTracks));
      if (!t)
        continue;
      if (t->Pt()<fMinJetTrackPt) 
        continue;

      Int_t index = 1;
      if(fMC && t->Charge() == 0)
	index =- 1;

      fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), (iTracks + 123) * index);
    }
  }

  if (clus) {
    Double_t vertex[3] = {0, 0, 0};
    fEvent->GetPrimaryVertex()->GetXYZ(vertex);
    const Int_t Nclus = clus->GetEntries();
    for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
      AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
      if (!c)
        continue;
      if (!c->IsEMCAL())
        continue;
      TLorentzVector nPart;
      c->GetMomentum(nPart, vertex);
      Double_t et = nPart.Pt();
      if (et < fMinJetClusPt) 
        continue;
      fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 123);
    }
  }

  // run jet finder
  fjw.Run();

  // get geometry
  AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
  if (!geom) {
    AliFatal("Can not create geometry");
    return;
  }
  
  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());
    Double_t vertex[3] = {0, 0, 0};
    fEvent->GetPrimaryVertex()->GetXYZ(vertex);
    vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
    Int_t nt            = 0;
    Int_t nc            = 0;
    Double_t neutralE   = 0;
    Double_t maxTrack   = 0;
    Double_t maxCluster = 0;
    Int_t gall          = 0;
    Int_t gemc          = 0;

    jet->SetNumberOfTracks(constituents.size());
    jet->SetNumberOfClusters(constituents.size());

    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() * TMath::RadToDeg();
        Double_t geta = constituents[ic].eta();
        if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
            (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
          ++gemc;
        continue;
      }	else if (fMC || uid >= 0) {
	Int_t tid = TMath::Abs(uid) - 123;
        AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
        if (t) {
	  if (uid >= 0)
	    neutralE += t->P();
          if (t->Pt() > maxTrack)
            maxTrack = t->Pt();
          jet->AddTrackAt(tid, nt);
          nt++;
        }
      } else {
	Int_t cid = TMath::Abs(uid) - 123;
	AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
	if (c) {
	  TLorentzVector nP;
	  c->GetMomentum(nP, vertex);
	  neutralE += nP.P();
	  if (nP.Pt() > maxCluster)
	    maxCluster = nP.Pt();
	  jet->AddClusterAt(cid, nc);
	  nc++;
	}
      }
    }
    jet->SetNumberOfTracks(nt);
    jet->SetNumberOfClusters(nc);
    jet->SortConstituents();
    jet->SetMaxTrackPt(maxTrack);
    jet->SetMaxClusterPt(maxCluster);
    jet->SetNEF(neutralE / jet->E());
    jet->SetArea(fjw.GetJetArea(ij));
    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++;
  }
}
Commit	Line	Data
48d874ff	1	// $Id$
542ac102	2	//
	3	// Emcal jet finder task.
	4	//
b65fac7a	5	// Authors: C.Loizides, S.Aiola
48d874ff	6
6ea168d0	7	#include <TChain.h>
48d874ff	8	#include <TClonesArray.h>
6ea168d0	9	#include <TList.h>
6ea168d0	10	#include <TLorentzVector.h>
b65fac7a	11
48d874ff	12	#include "AliAnalysisManager.h"
6ea168d0	13	#include "AliCentrality.h"
f5f3c8e9	14	#include "AliEMCALGeometry.h"
f5f3c8e9	15	#include "AliEmcalJet.h"
63206072	16	#include "AliFJWrapper.h"
35789a2d	17	#include "AliVCluster.h"
b65fac7a	18	#include "AliVParticle.h"
	19	#include "AliVEvent.h"
	20	#include "AliMCEvent.h"
	21
	22	#include "AliEmcalJetTask.h"
48d874ff	23
914d486c	24	ClassImp(AliEmcalJetTask)
48d874ff	25
48d874ff	26	//________________________________________________________________________
914d486c	27	AliEmcalJetTask::AliEmcalJetTask() :
914d486c	28	AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0	29	fTracksName("Tracks"),
	30	fCaloName("CaloClusters"),
	31	fJetsName("Jets"),
b65fac7a	32	fMC(kFALSE),
6ea168d0	33	fAlgo(1),
	34	fRadius(0.4),
	35	fType(0),
914d486c	36	fMinJetTrackPt(0.15),
914d486c	37	fMinJetClusPt(0.15),
f5f3c8e9	38	fMinJetArea(0.01),
f5f3c8e9	39	fMinJetPt(1.0),
b65fac7a	40	fJets(0),
b65fac7a	41	fEvent(0)
6ea168d0	42	{
6ea168d0	43	// Default constructor.
6ea168d0	44	}
	45
	46	//________________________________________________________________________
914d486c	47	AliEmcalJetTask::AliEmcalJetTask(const char *name) :
63206072	48	AliAnalysisTaskSE(name),
7efbea04	49	fTracksName("Tracks"),
48d874ff	50	fCaloName("CaloClusters"),
7efbea04	51	fJetsName("Jets"),
b65fac7a	52	fMC(kFALSE),
48d874ff	53	fAlgo(1),
	54	fRadius(0.4),
	55	fType(0),
6ea168d0	56	fMinJetTrackPt(0.15),
914d486c	57	fMinJetClusPt(0.15),
f5f3c8e9	58	fMinJetArea(0.01),
f5f3c8e9	59	fMinJetPt(1.0),
b65fac7a	60	fJets(0),
b65fac7a	61	fEvent(0)
48d874ff	62	{
48d874ff	63	// Standard constructor.
7efbea04	64
7efbea04	65	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff	66	}
	67
	68	//________________________________________________________________________
914d486c	69	AliEmcalJetTask::~AliEmcalJetTask()
48d874ff	70	{
	71	// Destructor
	72	}
	73
	74	//________________________________________________________________________
914d486c	75	void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff	76	{
	77	// Create user objects.
	78
914d486c	79	fJets = new TClonesArray("AliEmcalJet");
48d874ff	80	fJets->SetName(fJetsName);
48d874ff	81	}
	82
	83	//________________________________________________________________________
914d486c	84	void AliEmcalJetTask::UserExec(Option_t *)
48d874ff	85	{
	86	// Main loop, called for each event.
	87
b65fac7a	88	// get the event
	89	if (fMC)
	90	fEvent = MCEvent();
	91	else
	92	fEvent = InputEvent();
	93
63206072	94	// add jets to event if not yet there
b65fac7a	95	if (!(fEvent->FindListObject(fJetsName)))
b65fac7a	96	fEvent->AddObject(fJets);
48d874ff	97
63206072	98	// delete jet output
	99	fJets->Delete();
	100
b65fac7a	101	if (!fEvent) {
	102	AliError(Form("Could not get the event! fMC = %d", fMC));
	103	return;
	104	}
	105
63206072	106	// get input collections
48d874ff	107	TClonesArray *tracks = 0;
48d874ff	108	TClonesArray *clus = 0;
63206072	109	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
6ea168d0	110
b65fac7a	111	if ((fType == 0) \|\| (fType == 1)) {
7efbea04	112	if (fTracksName == "Tracks")
48d874ff	113	am->LoadBranch("Tracks");
b65fac7a	114	tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
48d874ff	115	if (!tracks) {
653159b9	116	AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
48d874ff	117	return;
	118	}
	119	}
b65fac7a	120	if ((fType == 0) \|\| (fType == 2)) {
48d874ff	121	if (fCaloName == "CaloClusters")
48d874ff	122	am->LoadBranch("CaloClusters");
b65fac7a	123	clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
48d874ff	124	if (!clus) {
653159b9	125	AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
48d874ff	126	return;
48d874ff	127	}
48d874ff	128	}
63206072	129
	130	// get centrality
	131	Float_t cent = -1;
b65fac7a	132	AliCentrality *centrality = fEvent->GetCentrality();
63206072	133	if (centrality)
	134	cent = centrality->GetCentralityPercentile("V0M");
	135	else
	136	cent=99; // probably pp data
b65fac7a	137	if (cent < 0) {
63206072	138	AliError(Form("Centrality negative: %f", cent));
	139	return;
	140	}
653159b9	141
6ea168d0	142	FindJets(tracks, clus, fAlgo, fRadius, cent);
48d874ff	143	}
	144
	145	//________________________________________________________________________
914d486c	146	void AliEmcalJetTask::Terminate(Option_t *)
48d874ff	147	{
48d874ff	148	// Called once at the end of the analysis.
48d874ff	149	}
	150
	151	//________________________________________________________________________
914d486c	152	void AliEmcalJetTask::FindJets(TObjArray tracks, TObjArray clus, Int_t algo, Double_t radius, Float_t /cent/)
48d874ff	153	{
	154	// Find jets.
	155
d03084cd	156	if (!tracks && !clus)
	157	return;
	158
48d874ff	159	TString name("kt");
	160	fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
	161	if (algo>=1) {
	162	name = "antikt";
	163	jalgo = fastjet::antikt_algorithm;
	164	}
	165
	166	AliFJWrapper fjw(name, name);
f5f3c8e9	167	fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
48d874ff	168	fjw.SetR(radius);
6ea168d0	169	fjw.SetAlgorithm(jalgo);
48d874ff	170	fjw.SetMaxRap(0.9);
	171	fjw.Clear();
	172
	173	if (tracks) {
	174	const Int_t Ntracks = tracks->GetEntries();
	175	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
b65fac7a	176	AliVParticle t = static_cast<AliVParticle>(tracks->At(iTracks));
48d874ff	177	if (!t)
48d874ff	178	continue;
914d486c	179	if (t->Pt()<fMinJetTrackPt)
914d486c	180	continue;
b65fac7a	181
	182	Int_t index = 1;
	183	if(fMC && t->Charge() == 0)
	184	index =- 1;
	185
	186	fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), (iTracks + 123) * index);
48d874ff	187	}
48d874ff	188	}
6ea168d0	189
48d874ff	190	if (clus) {
48d874ff	191	Double_t vertex[3] = {0, 0, 0};
b65fac7a	192	fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff	193	const Int_t Nclus = clus->GetEntries();
914d486c	194	for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
7efbea04	195	AliVCluster c = dynamic_cast<AliVCluster>(clus->At(iClus));
d03084cd	196	if (!c)
d03084cd	197	continue;
48d874ff	198	if (!c->IsEMCAL())
	199	continue;
	200	TLorentzVector nPart;
	201	c->GetMomentum(nPart, vertex);
63206072	202	Double_t et = nPart.Pt();
b65fac7a	203	if (et < fMinJetClusPt)
914d486c	204	continue;
b65fac7a	205	fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 123);
48d874ff	206	}
	207	}
	208
7efbea04	209	// run jet finder
48d874ff	210	fjw.Run();
f5f3c8e9	211
	212	// get geometry
	213	AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
	214	if (!geom) {
	215	AliFatal("Can not create geometry");
	216	return;
	217	}
6ea168d0	218
48d874ff	219	std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
b65fac7a	220	for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9	221	if (jets_incl[ij].perp()<fMinJetPt)
	222	continue;
	223	if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff	224	continue;
914d486c	225	AliEmcalJet jet = new ((fJets)[jetCount])
914d486c	226	AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
6ea168d0	227	Double_t vertex[3] = {0, 0, 0};
b65fac7a	228	fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff	229	vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
b65fac7a	230	Int_t nt = 0;
b65fac7a	231	Int_t nc = 0;
f5f3c8e9	232	Double_t neutralE = 0;
	233	Double_t maxTrack = 0;
	234	Double_t maxCluster = 0;
	235	Int_t gall = 0;
	236	Int_t gemc = 0;
b65fac7a	237
f5f3c8e9	238	jet->SetNumberOfTracks(constituents.size());
f5f3c8e9	239	jet->SetNumberOfClusters(constituents.size());
b65fac7a	240
b65fac7a	241	for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff	242	Int_t uid = constituents[ic].user_index();
b65fac7a	243
b65fac7a	244	if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
f5f3c8e9	245	++gall;
b65fac7a	246	Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
f5f3c8e9	247	Double_t geta = constituents[ic].eta();
b65fac7a	248	if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
b65fac7a	249	(geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
f5f3c8e9	250	++gemc;
f5f3c8e9	251	continue;
b65fac7a	252	} else if (fMC \|\| uid >= 0) {
	253	Int_t tid = TMath::Abs(uid) - 123;
	254	AliVParticle t = static_cast<AliVParticle>(tracks->At(tid));
d03084cd	255	if (t) {
b65fac7a	256	if (uid >= 0)
	257	neutralE += t->P();
	258	if (t->Pt() > maxTrack)
	259	maxTrack = t->Pt();
	260	jet->AddTrackAt(tid, nt);
d03084cd	261	nt++;
d03084cd	262	}
6ea168d0	263	} else {
b65fac7a	264	Int_t cid = TMath::Abs(uid) - 123;
	265	AliVCluster c = static_cast<AliVCluster>(clus->At(cid));
	266	if (c) {
	267	TLorentzVector nP;
	268	c->GetMomentum(nP, vertex);
	269	neutralE += nP.P();
	270	if (nP.Pt() > maxCluster)
	271	maxCluster = nP.Pt();
	272	jet->AddClusterAt(cid, nc);
	273	nc++;
	274	}
48d874ff	275	}
48d874ff	276	}
35789a2d	277	jet->SetNumberOfTracks(nt);
35789a2d	278	jet->SetNumberOfClusters(nc);
f5f3c8e9	279	jet->SortConstituents();
6ea168d0	280	jet->SetMaxTrackPt(maxTrack);
6ea168d0	281	jet->SetMaxClusterPt(maxCluster);
b65fac7a	282	jet->SetNEF(neutralE / jet->E());
f5f3c8e9	283	jet->SetArea(fjw.GetJetArea(ij));
b65fac7a	284	if (gall > 0)
b65fac7a	285	jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9	286	else
f5f3c8e9	287	jet->SetAreaEmc(-1);
b65fac7a	288	if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
	289	(jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
	290	(jet->Eta() > geom->GetArm1EtaMin()) &&
	291	(jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9	292	jet->SetAxisInEmcal(kTRUE);
48d874ff	293	jetCount++;
	294	}
	295	}