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