[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliEmcalIsolatedPhotonsTask.cxx

// $Id: AliEmcalIsolatedPhotonsTask.cxx  $

#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TParticle.h>

#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliVCluster.h"
#include "AliESDtrack.h"
#include "AliEmcalJet.h"
#include "AliAODTrack.h"
#include "AliEmcalJet.h"
#include "AliVEventHandler.h"
#include "AliPicoTrack.h"

#include "AliEmcalIsolatedPhotonsTask.h"

ClassImp(AliEmcalIsolatedPhotonsTask)

//________________________________________________________________________
AliEmcalIsolatedPhotonsTask::AliEmcalIsolatedPhotonsTask() : 
  AliAnalysisTaskSE("AliEmcalIsolatedPhotonsTask"),
  fOutput(0),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fTrgClusName("ClustersL1GAMMAFEE"),
  fTracks(0),
  fCaloClusters(0),
  fJets(0),
  fTrgClusters(0),
  fCent(0),
  fHistCentrality(0),
  Ptbins(400),
  Ptlow(0),
  Ptup(200),
  Ebins(400),
  Elow(0),
  Eup(200)
{
  // Default constructor.
  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";

  for (Int_t i = 0; i < 4; i++) {
    fHistJetsE[i] = 0;
    fHistJetsNE[i] = 0;
    fHistJetsNEF[i] = 0;
    fHistJetsZ[i] = 0;
    fHistLeadingJetE[i] = 0;
    fHistTracksPtLJ[i] = 0;
    fHistClusELJ[i] = 0;
    fHistTracksPtBkg[i] = 0;
    fHistClusEBkg[i] = 0;
  }
}

//________________________________________________________________________
AliEmcalIsolatedPhotonsTask::AliEmcalIsolatedPhotonsTask(const char *name) : 
  AliAnalysisTaskSE("AliEmcalIsolatedPhotonsTask"),
  fOutput(0),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fTrgClusName("ClustersL1GAMMAFEE"),
  fTracks(0),
  fCaloClusters(0),
  fJets(0),
  fTrgClusters(0),
  fCent(0),
  fHistCentrality(0),
  Ptbins(400),
  Ptlow(0),
  Ptup(200),
  Ebins(400),
  Elow(0),
  Eup(200)
{
  // Standard constructor.

  if (!name)
    return;

  SetName(name);
  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";

  for (Int_t i = 0; i < 4; i++) {
    fHistJetsE[i] = 0;
    fHistJetsNE[i] = 0;
    fHistJetsNEF[i] = 0;
    fHistJetsZ[i] = 0;
    fHistLeadingJetE[i] = 0;
    fHistTracksPtLJ[i] = 0;
    fHistClusELJ[i] = 0;
    fHistTracksPtBkg[i] = 0;
    fHistClusEBkg[i] = 0;
  }

  DefineInput(0,TChain::Class());
  DefineOutput(1,TList::Class());
}

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

//________________________________________________________________________
void AliEmcalIsolatedPhotonsTask::UserCreateOutputObjects()
{
  // Create histograms
  
  fOutput = new TList();
  fOutput->SetOwner();  // IMPORTANT!
  
  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", Ebins, 0, 100);
  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
  fHistCentrality->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistCentrality);

  TString histname;

  for (Int_t i = 0; i < 4; i++) {
    histname = "fHistJetsE_";
    histname += i;
    fHistJetsE[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
    fHistJetsE[i]->GetXaxis()->SetTitle("E [GeV]");
    fHistJetsE[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistJetsE[i]);
    
    histname = "fHistJetsNE_";
    histname += i;
    fHistJetsNE[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
    fHistJetsNE[i]->GetXaxis()->SetTitle("E [GeV]");
    fHistJetsNE[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistJetsNE[i]);
    
    histname = "fHistJetsNEF_";
    histname += i;
    fHistJetsNEF[i] = new TH1F(histname.Data(), histname.Data(), Ebins, 0, 1.2);
    fHistJetsNEF[i]->GetXaxis()->SetTitle("NEF");
    fHistJetsNEF[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistJetsNEF[i]);

    histname = "fHistJetsZ_";
    histname += i;
    fHistJetsZ[i] = new TH1F(histname.Data(), histname.Data(), Ebins, 0, 1.2);
    fHistJetsZ[i]->GetXaxis()->SetTitle("Z");
    fHistJetsZ[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistJetsZ[i]);

    histname = "fHistLeadingJetE_";
    histname += i;
    fHistLeadingJetE[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
    fHistLeadingJetE[i]->GetXaxis()->SetTitle("E [GeV]");
    fHistLeadingJetE[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistLeadingJetE[i]);
    
    histname = "fHistTracksPtLJ_";
    histname += i;
    fHistTracksPtLJ[i] = new TH1F(histname.Data(), histname.Data(), Ptbins, Ptlow, Ptup);
    fHistTracksPtLJ[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
    fHistTracksPtLJ[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistTracksPtLJ[i]);
    
    histname = "fHistClusELJ_";
    histname += i;
    fHistClusELJ[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
    fHistClusELJ[i]->GetXaxis()->SetTitle("E [GeV]");
    fHistClusELJ[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistClusELJ[i]);

    histname = "fHistTracksPtBkg_";
    histname += i;
    fHistTracksPtBkg[i] = new TH1F(histname.Data(), histname.Data(), Ptbins, Ptlow, Ptup);
    fHistTracksPtBkg[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
    fHistTracksPtBkg[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistTracksPtBkg[i]);
    
    histname = "fHistClusEBkg_";
    histname += i;
    fHistClusEBkg[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
    fHistClusEBkg[i]->GetXaxis()->SetTitle("E [GeV]");
    fHistClusEBkg[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistClusEBkg[i]);
  }

  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
}

void AliEmcalIsolatedPhotonsTask::RetrieveEventObjects()
{
  fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
  if (!fCaloClusters) {
    AliWarning(Form("Could not retrieve clusters!")); 
  }

  fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
  if (!fTracks) {
    AliWarning(Form("Could not retrieve tracks!")); 
  }

  fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
  if (!fJets) {
    AliWarning(Form("Could not retrieve jets!")); 
  }

  if (strcmp(fTrgClusName,"")) {
    fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
    if (!fTrgClusters) {
      AliWarning(Form("Could not retrieve trigger clusters!")); 
    }
  }

  fCent = InputEvent()->GetCentrality();
}

AliVTrack* AliEmcalIsolatedPhotonsTask::GetTrack(const Int_t i) const
{
  if (fTracks)
    return dynamic_cast<AliVTrack*>(fTracks->At(i));
  else
    return 0;
}

Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfTracks() const
{
  if (fTracks)
    return fTracks->GetEntriesFast();
  else
    return 0;
}

AliVCluster* AliEmcalIsolatedPhotonsTask::GetCaloCluster(const Int_t i) const
{ 
  if (fCaloClusters)
    return dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
  else
    return 0;
}

Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfCaloClusters() const
{ 
  if (fCaloClusters)
    return fCaloClusters->GetEntriesFast();
  else
    return 0;
}

AliEmcalJet* AliEmcalIsolatedPhotonsTask::GetJet(const Int_t i) const
{
  if (fJets)
    return dynamic_cast<AliEmcalJet*>(fJets->At(i));
  else
    return 0;
}

Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfJets() const
{
  if (fJets)
    return fJets->GetEntriesFast();
  else
    return 0;
}

AliVCluster* AliEmcalIsolatedPhotonsTask::GetTrgCluster(const Int_t i) const
{
  if (fTrgClusters)
    return dynamic_cast<AliVCluster*>(fTrgClusters->At(i));
  else
    return 0;
}

Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfTrgClusters() const
{
  if (fTrgClusters)
    return fTrgClusters->GetEntriesFast();
  else
    return 0;
}

void AliEmcalIsolatedPhotonsTask::FillHistograms()
{
  Float_t cent = fCent->GetCentralityPercentile("V0M");

  fHistCentrality->Fill(cent);

  Int_t centbin=-1;
  if(cent>=0 && cent<10) centbin=0;
  else if(cent>=10 && cent<30) centbin=1;
  else if(cent>=30 && cent<50) centbin=2;
  else if(cent>=50 && cent<=100) centbin=3;

  // Jet loop
  Int_t njets =  GetNumberOfJets();
  //cout << njets << " jets" << endl;
  Float_t maxJetEnergy = 0;
  Int_t maxJetIndex = -1;
  for (Int_t ij = 0; ij < njets; ij++) {

    AliEmcalJet* jet = GetJet(ij);

    if (!jet) {
      printf("ERROR: Could not receive jet %d\n", ij);
      continue;
    }  
    
    if (jet->E() <= 0)
      continue;

    fHistJetsE[centbin]->Fill(jet->E());
    fHistJetsNEF[centbin]->Fill(jet->NEF());
    fHistJetsNE[centbin]->Fill(jet->E() * jet->NEF());

    //cout << "****** jet id = " << ij << " ******" << endl;

    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
      Int_t trackid = jet->TrackAt(it);
      AliVTrack *track = GetTrack(trackid);
      //cout << "track id = " << trackid << endl;
      if (track)
	fHistJetsZ[centbin]->Fill(track->Pt() / jet->E());
    }
    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
      Int_t clusterid = jet->ClusterAt(ic);
      //cout << "cluster id = " << clusterid << endl;
      AliVCluster *cluster = GetCaloCluster(clusterid);
      if (cluster)
	fHistJetsZ[centbin]->Fill(cluster->E() / jet->E());
    }

    if (maxJetEnergy < jet->E()) {
      maxJetEnergy = jet->E();
      maxJetIndex = ij;
    }
  } //jet loop 

  
  if (!(maxJetEnergy > 0) || maxJetIndex < 0) 
    return;

  fHistLeadingJetE[centbin]->Fill(maxJetEnergy);
  
  AliEmcalJet* jet = GetJet(maxJetIndex);
  if (!jet)
    return;
  
  jet->SortConstituents();
  
  // Cluster loop
  Int_t clusJetId = 0;
  Int_t nclusters =  GetNumberOfCaloClusters();
  for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
    AliVCluster* cluster = GetCaloCluster(iClusters);
    if (!cluster) {
      printf("ERROR: Could not receive cluster %d\n", iClusters);
      continue;
    }  
    
    if (!(cluster->IsEMCAL())) continue;
    
    if (jet->GetNumberOfClusters() > 0) {
      if (jet->ClusterAt(clusJetId) < iClusters) {
	if (clusJetId < jet->GetNumberOfClusters() - 1)
	  clusJetId++;
      }
      
      if (jet->ClusterAt(clusJetId) == iClusters) {
	fHistClusELJ[centbin]->Fill(cluster->E());
      }
      else {
	fHistClusEBkg[centbin]->Fill(cluster->E());
      }
    }
    else {
      fHistClusEBkg[centbin]->Fill(cluster->E());
    }
  } //cluster loop 
  
  
  // Track loop 
  Int_t trackJetId = 0;
  Int_t ntracks = GetNumberOfTracks();
  for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
    AliVTrack* track = GetTrack(iTracks);         
    if(!track) {
      AliError(Form("ERROR: Could not retrieve track %d",iTracks)); 
      continue; 
    }
    
    if (!AcceptTrack(track)) continue;
    
    if (jet->GetNumberOfTracks() > 0) {
      if (jet->TrackAt(trackJetId) < iTracks) {
	if (trackJetId < jet->GetNumberOfTracks() - 1)
	  trackJetId++;
      }
      
      if (jet->TrackAt(trackJetId) == iTracks) {
	fHistTracksPtLJ[centbin]->Fill(track->Pt());
      }
      else {
	fHistTracksPtBkg[centbin]->Fill(track->Pt());
      }
    }
    else {
      fHistTracksPtBkg[centbin]->Fill(track->Pt());
    }
  }
}


//________________________________________________________________________
Bool_t AliEmcalIsolatedPhotonsTask::AcceptTrack(AliVTrack* /*track*/)
{
  return kTRUE;
}

//________________________________________________________________________
void AliEmcalIsolatedPhotonsTask::UserExec(Option_t *) 
{
  // Main loop, called for each event.
  // Add jets to event if not yet there

  RetrieveEventObjects();

  FillHistograms();
    
  // information for this iteration of the UserExec in the container
  PostData(1, fOutput);
}

//________________________________________________________________________
void AliEmcalIsolatedPhotonsTask::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
25283b37	1	// $Id: AliEmcalIsolatedPhotonsTask.cxx $
	2
	3	#include <TChain.h>
	4	#include <TClonesArray.h>
	5	#include <TH1F.h>
	6	#include <TH2F.h>
	7	#include <TList.h>
	8	#include <TLorentzVector.h>
	9	#include <TParticle.h>
	10
	11	#include "AliAnalysisManager.h"
	12	#include "AliCentrality.h"
f0a0fd33	13	#include "AliVCluster.h"
25283b37	14	#include "AliESDtrack.h"
25283b37	15	#include "AliEmcalJet.h"
25283b37	16	#include "AliAODTrack.h"
25283b37	17	#include "AliEmcalJet.h"
25283b37	18	#include "AliVEventHandler.h"
35789a2d	19	#include "AliPicoTrack.h"
25283b37	20
	21	#include "AliEmcalIsolatedPhotonsTask.h"
	22
	23	ClassImp(AliEmcalIsolatedPhotonsTask)
	24
	25	//________________________________________________________________________
	26	AliEmcalIsolatedPhotonsTask::AliEmcalIsolatedPhotonsTask() :
	27	AliAnalysisTaskSE("AliEmcalIsolatedPhotonsTask"),
	28	fOutput(0),
	29	fTracksName("Tracks"),
	30	fCaloName("CaloClusters"),
	31	fJetsName("Jets"),
e82e282c	32	fTrgClusName("ClustersL1GAMMAFEE"),
25283b37	33	fTracks(0),
	34	fCaloClusters(0),
	35	fJets(0),
e82e282c	36	fTrgClusters(0),
f0a0fd33	37	fCent(0),
	38	fHistCentrality(0),
	39	Ptbins(400),
25283b37	40	Ptlow(0),
f0a0fd33	41	Ptup(200),
f0a0fd33	42	Ebins(400),
25283b37	43	Elow(0),
f0a0fd33	44	Eup(200)
25283b37	45	{
	46	// Default constructor.
	47	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
f0a0fd33	48
	49	for (Int_t i = 0; i < 4; i++) {
	50	fHistJetsE[i] = 0;
	51	fHistJetsNE[i] = 0;
	52	fHistJetsNEF[i] = 0;
	53	fHistJetsZ[i] = 0;
	54	fHistLeadingJetE[i] = 0;
	55	fHistTracksPtLJ[i] = 0;
	56	fHistClusELJ[i] = 0;
	57	fHistTracksPtBkg[i] = 0;
	58	fHistClusEBkg[i] = 0;
	59	}
25283b37	60	}
	61
	62	//________________________________________________________________________
	63	AliEmcalIsolatedPhotonsTask::AliEmcalIsolatedPhotonsTask(const char *name) :
	64	AliAnalysisTaskSE("AliEmcalIsolatedPhotonsTask"),
	65	fOutput(0),
	66	fTracksName("Tracks"),
	67	fCaloName("CaloClusters"),
	68	fJetsName("Jets"),
e82e282c	69	fTrgClusName("ClustersL1GAMMAFEE"),
25283b37	70	fTracks(0),
	71	fCaloClusters(0),
	72	fJets(0),
e82e282c	73	fTrgClusters(0),
f0a0fd33	74	fCent(0),
	75	fHistCentrality(0),
	76	Ptbins(400),
25283b37	77	Ptlow(0),
f0a0fd33	78	Ptup(200),
f0a0fd33	79	Ebins(400),
25283b37	80	Elow(0),
f0a0fd33	81	Eup(200)
25283b37	82	{
	83	// Standard constructor.
	84
	85	if (!name)
	86	return;
	87
	88	SetName(name);
	89	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
	90
f0a0fd33	91	for (Int_t i = 0; i < 4; i++) {
	92	fHistJetsE[i] = 0;
	93	fHistJetsNE[i] = 0;
	94	fHistJetsNEF[i] = 0;
	95	fHistJetsZ[i] = 0;
	96	fHistLeadingJetE[i] = 0;
	97	fHistTracksPtLJ[i] = 0;
	98	fHistClusELJ[i] = 0;
	99	fHistTracksPtBkg[i] = 0;
	100	fHistClusEBkg[i] = 0;
	101	}
	102
25283b37	103	DefineInput(0,TChain::Class());
	104	DefineOutput(1,TList::Class());
	105	}
	106
	107	//________________________________________________________________________
	108	AliEmcalIsolatedPhotonsTask::~AliEmcalIsolatedPhotonsTask()
	109	{
	110	// Destructor
	111	}
	112
	113	//________________________________________________________________________
	114	void AliEmcalIsolatedPhotonsTask::UserCreateOutputObjects()
	115	{
f0a0fd33	116	// Create histograms
25283b37	117
	118	fOutput = new TList();
	119	fOutput->SetOwner(); // IMPORTANT!
	120
f0a0fd33	121	fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", Ebins, 0, 100);
	122	fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
	123	fHistCentrality->GetYaxis()->SetTitle("counts");
	124	fOutput->Add(fHistCentrality);
25283b37	125
f0a0fd33	126	TString histname;
	127
	128	for (Int_t i = 0; i < 4; i++) {
	129	histname = "fHistJetsE_";
	130	histname += i;
	131	fHistJetsE[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
	132	fHistJetsE[i]->GetXaxis()->SetTitle("E [GeV]");
	133	fHistJetsE[i]->GetYaxis()->SetTitle("counts");
	134	fOutput->Add(fHistJetsE[i]);
25283b37	135
f0a0fd33	136	histname = "fHistJetsNE_";
	137	histname += i;
	138	fHistJetsNE[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
	139	fHistJetsNE[i]->GetXaxis()->SetTitle("E [GeV]");
	140	fHistJetsNE[i]->GetYaxis()->SetTitle("counts");
	141	fOutput->Add(fHistJetsNE[i]);
25283b37	142
f0a0fd33	143	histname = "fHistJetsNEF_";
	144	histname += i;
	145	fHistJetsNEF[i] = new TH1F(histname.Data(), histname.Data(), Ebins, 0, 1.2);
	146	fHistJetsNEF[i]->GetXaxis()->SetTitle("NEF");
	147	fHistJetsNEF[i]->GetYaxis()->SetTitle("counts");
	148	fOutput->Add(fHistJetsNEF[i]);
	149
	150	histname = "fHistJetsZ_";
	151	histname += i;
	152	fHistJetsZ[i] = new TH1F(histname.Data(), histname.Data(), Ebins, 0, 1.2);
	153	fHistJetsZ[i]->GetXaxis()->SetTitle("Z");
	154	fHistJetsZ[i]->GetYaxis()->SetTitle("counts");
	155	fOutput->Add(fHistJetsZ[i]);
	156
	157	histname = "fHistLeadingJetE_";
	158	histname += i;
	159	fHistLeadingJetE[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
	160	fHistLeadingJetE[i]->GetXaxis()->SetTitle("E [GeV]");
	161	fHistLeadingJetE[i]->GetYaxis()->SetTitle("counts");
	162	fOutput->Add(fHistLeadingJetE[i]);
25283b37	163
f0a0fd33	164	histname = "fHistTracksPtLJ_";
	165	histname += i;
	166	fHistTracksPtLJ[i] = new TH1F(histname.Data(), histname.Data(), Ptbins, Ptlow, Ptup);
	167	fHistTracksPtLJ[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
	168	fHistTracksPtLJ[i]->GetYaxis()->SetTitle("counts");
	169	fOutput->Add(fHistTracksPtLJ[i]);
	170
	171	histname = "fHistClusELJ_";
	172	histname += i;
	173	fHistClusELJ[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
	174	fHistClusELJ[i]->GetXaxis()->SetTitle("E [GeV]");
	175	fHistClusELJ[i]->GetYaxis()->SetTitle("counts");
	176	fOutput->Add(fHistClusELJ[i]);
	177
	178	histname = "fHistTracksPtBkg_";
	179	histname += i;
	180	fHistTracksPtBkg[i] = new TH1F(histname.Data(), histname.Data(), Ptbins, Ptlow, Ptup);
	181	fHistTracksPtBkg[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
	182	fHistTracksPtBkg[i]->GetYaxis()->SetTitle("counts");
	183	fOutput->Add(fHistTracksPtBkg[i]);
	184
	185	histname = "fHistClusEBkg_";
	186	histname += i;
	187	fHistClusEBkg[i] = new TH1F(histname.Data(), histname.Data(), Ebins, Elow, Eup);
	188	fHistClusEBkg[i]->GetXaxis()->SetTitle("E [GeV]");
	189	fHistClusEBkg[i]->GetYaxis()->SetTitle("counts");
	190	fOutput->Add(fHistClusEBkg[i]);
e82e282c	191	}
e82e282c	192
25283b37	193	PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
	194	}
	195
	196	void AliEmcalIsolatedPhotonsTask::RetrieveEventObjects()
	197	{
	198	fCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
f0a0fd33	199	if (!fCaloClusters) {
	200	AliWarning(Form("Could not retrieve clusters!"));
	201	}
25283b37	202
f0a0fd33	203	fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
25283b37	204	if (!fTracks) {
e82e282c	205	AliWarning(Form("Could not retrieve tracks!"));
25283b37	206	}
f0a0fd33	207
f0a0fd33	208	fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
25283b37	209	if (!fJets) {
e82e282c	210	AliWarning(Form("Could not retrieve jets!"));
e82e282c	211	}
f0a0fd33	212
	213	if (strcmp(fTrgClusName,"")) {
	214	fTrgClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
	215	if (!fTrgClusters) {
	216	AliWarning(Form("Could not retrieve trigger clusters!"));
	217	}
25283b37	218	}
f0a0fd33	219
f0a0fd33	220	fCent = InputEvent()->GetCentrality();
25283b37	221	}
	222
	223	AliVTrack* AliEmcalIsolatedPhotonsTask::GetTrack(const Int_t i) const
	224	{
e82e282c	225	if (fTracks)
	226	return dynamic_cast<AliVTrack*>(fTracks->At(i));
	227	else
	228	return 0;
25283b37	229	}
	230
	231	Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfTracks() const
	232	{
e82e282c	233	if (fTracks)
	234	return fTracks->GetEntriesFast();
	235	else
	236	return 0;
25283b37	237	}
	238
	239	AliVCluster* AliEmcalIsolatedPhotonsTask::GetCaloCluster(const Int_t i) const
e82e282c	240	{
	241	if (fCaloClusters)
	242	return dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
	243	else
	244	return 0;
25283b37	245	}
	246
	247	Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfCaloClusters() const
e82e282c	248	{
	249	if (fCaloClusters)
	250	return fCaloClusters->GetEntriesFast();
	251	else
	252	return 0;
25283b37	253	}
	254
	255	AliEmcalJet* AliEmcalIsolatedPhotonsTask::GetJet(const Int_t i) const
	256	{
e82e282c	257	if (fJets)
	258	return dynamic_cast<AliEmcalJet*>(fJets->At(i));
	259	else
	260	return 0;
25283b37	261	}
	262
	263	Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfJets() const
	264	{
e82e282c	265	if (fJets)
	266	return fJets->GetEntriesFast();
	267	else
	268	return 0;
	269	}
	270
	271	AliVCluster* AliEmcalIsolatedPhotonsTask::GetTrgCluster(const Int_t i) const
	272	{
	273	if (fTrgClusters)
	274	return dynamic_cast<AliVCluster*>(fTrgClusters->At(i));
	275	else
	276	return 0;
	277	}
	278
	279	Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfTrgClusters() const
	280	{
	281	if (fTrgClusters)
	282	return fTrgClusters->GetEntriesFast();
	283	else
	284	return 0;
25283b37	285	}
	286
	287	void AliEmcalIsolatedPhotonsTask::FillHistograms()
	288	{
f0a0fd33	289	Float_t cent = fCent->GetCentralityPercentile("V0M");
25283b37	290
f0a0fd33	291	fHistCentrality->Fill(cent);
e82e282c	292
f0a0fd33	293	Int_t centbin=-1;
	294	if(cent>=0 && cent<10) centbin=0;
	295	else if(cent>=10 && cent<30) centbin=1;
	296	else if(cent>=30 && cent<50) centbin=2;
	297	else if(cent>=50 && cent<=100) centbin=3;
25283b37	298
	299	// Jet loop
	300	Int_t njets = GetNumberOfJets();
	301	//cout << njets << " jets" << endl;
f0a0fd33	302	Float_t maxJetEnergy = 0;
f0a0fd33	303	Int_t maxJetIndex = -1;
25283b37	304	for (Int_t ij = 0; ij < njets; ij++) {
f0a0fd33	305
25283b37	306	AliEmcalJet* jet = GetJet(ij);
f0a0fd33	307
25283b37	308	if (!jet) {
	309	printf("ERROR: Could not receive jet %d\n", ij);
	310	continue;
	311	}
	312
f0a0fd33	313	if (jet->E() <= 0)
	314	continue;
	315
	316	fHistJetsE[centbin]->Fill(jet->E());
	317	fHistJetsNEF[centbin]->Fill(jet->NEF());
	318	fHistJetsNE[centbin]->Fill(jet->E() * jet->NEF());
	319
	320	//cout << "**** jet id = " << ij << " ****" << endl;
	321
35789a2d	322	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
	323	Int_t trackid = jet->TrackAt(it);
	324	AliVTrack *track = GetTrack(trackid);
f0a0fd33	325	//cout << "track id = " << trackid << endl;
35789a2d	326	if (track)
f0a0fd33	327	fHistJetsZ[centbin]->Fill(track->Pt() / jet->E());
35789a2d	328	}
	329	for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
	330	Int_t clusterid = jet->ClusterAt(ic);
f0a0fd33	331	//cout << "cluster id = " << clusterid << endl;
35789a2d	332	AliVCluster *cluster = GetCaloCluster(clusterid);
35789a2d	333	if (cluster)
f0a0fd33	334	fHistJetsZ[centbin]->Fill(cluster->E() / jet->E());
	335	}
	336
	337	if (maxJetEnergy < jet->E()) {
	338	maxJetEnergy = jet->E();
	339	maxJetIndex = ij;
35789a2d	340	}
25283b37	341	} //jet loop
e82e282c	342
f0a0fd33	343
	344	if (!(maxJetEnergy > 0) \|\| maxJetIndex < 0)
	345	return;
	346
	347	fHistLeadingJetE[centbin]->Fill(maxJetEnergy);
	348
	349	AliEmcalJet* jet = GetJet(maxJetIndex);
	350	if (!jet)
	351	return;
	352
	353	jet->SortConstituents();
	354
	355	// Cluster loop
	356	Int_t clusJetId = 0;
	357	Int_t nclusters = GetNumberOfCaloClusters();
	358	for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
	359	AliVCluster* cluster = GetCaloCluster(iClusters);
e82e282c	360	if (!cluster) {
	361	printf("ERROR: Could not receive cluster %d\n", iClusters);
	362	continue;
	363	}
	364
	365	if (!(cluster->IsEMCAL())) continue;
f0a0fd33	366
	367	if (jet->GetNumberOfClusters() > 0) {
	368	if (jet->ClusterAt(clusJetId) < iClusters) {
	369	if (clusJetId < jet->GetNumberOfClusters() - 1)
	370	clusJetId++;
	371	}
	372
	373	if (jet->ClusterAt(clusJetId) == iClusters) {
	374	fHistClusELJ[centbin]->Fill(cluster->E());
	375	}
	376	else {
	377	fHistClusEBkg[centbin]->Fill(cluster->E());
	378	}
	379	}
	380	else {
	381	fHistClusEBkg[centbin]->Fill(cluster->E());
	382	}
e82e282c	383	} //cluster loop
e82e282c	384
f0a0fd33	385
	386	// Track loop
	387	Int_t trackJetId = 0;
	388	Int_t ntracks = GetNumberOfTracks();
	389	for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
	390	AliVTrack* track = GetTrack(iTracks);
	391	if(!track) {
	392	AliError(Form("ERROR: Could not retrieve track %d",iTracks));
	393	continue;
25283b37	394	}
f0a0fd33	395
	396	if (!AcceptTrack(track)) continue;
	397
	398	if (jet->GetNumberOfTracks() > 0) {
	399	if (jet->TrackAt(trackJetId) < iTracks) {
	400	if (trackJetId < jet->GetNumberOfTracks() - 1)
	401	trackJetId++;
	402	}
25283b37	403
f0a0fd33	404	if (jet->TrackAt(trackJetId) == iTracks) {
	405	fHistTracksPtLJ[centbin]->Fill(track->Pt());
	406	}
	407	else {
	408	fHistTracksPtBkg[centbin]->Fill(track->Pt());
	409	}
	410	}
	411	else {
	412	fHistTracksPtBkg[centbin]->Fill(track->Pt());
25283b37	413	}
25283b37	414	}
f0a0fd33	415	}
	416
	417
	418	//________________________________________________________________________
	419	Bool_t AliEmcalIsolatedPhotonsTask::AcceptTrack(AliVTrack* /track/)
	420	{
	421	return kTRUE;
25283b37	422	}
	423
	424	//________________________________________________________________________
	425	void AliEmcalIsolatedPhotonsTask::UserExec(Option_t *)
	426	{
	427	// Main loop, called for each event.
	428	// Add jets to event if not yet there
	429
	430	RetrieveEventObjects();
	431
	432	FillHistograms();
	433
	434	// information for this iteration of the UserExec in the container
	435	PostData(1, fOutput);
	436	}
	437
	438	//________________________________________________________________________
	439	void AliEmcalIsolatedPhotonsTask::Terminate(Option_t *)
	440	{
	441	// Called once at the end of the analysis.
	442	}