[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 "AliESDCaloCluster.h"
#include "AliESDtrack.h"
#include "AliEmcalJet.h"
#include "AliFJWrapper.h"
#include "AliAODTrack.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.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"),
  fESDTrackCuts(0),
  fFilterBit(16),
  fTracks(0),
  fCaloClusters(0),
  fJets(0),
  fTrgClusters(0),
  fHistTracksPt(0),
  fHistClustersEnergy(0),
  fHistEPcorrelation(0),
  fHistJetsEnergy(0),
  fHistJetsNE(0),
  fHistJetsNEF(0),
  fHistJetsZ(0),
  fHistTrPhiEta(0),
  fHistClusPhiEta(0),
  fHistJetPhiEta(0),
  fHistMaxTrgCluster(0),
  Ptbins(100),
  Ptlow(0),
  Ptup(50),
  Ebins(100),
  Elow(0),
  Eup(50)
{
  // Default constructor.
  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
}

//________________________________________________________________________
AliEmcalIsolatedPhotonsTask::AliEmcalIsolatedPhotonsTask(const char *name) : 
  AliAnalysisTaskSE("AliEmcalIsolatedPhotonsTask"),
  fOutput(0),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fTrgClusName("ClustersL1GAMMAFEE"),
  fESDTrackCuts(0),
  fFilterBit(16),
  fTracks(0),
  fCaloClusters(0),
  fJets(0),
  fTrgClusters(0),
  fHistTracksPt(0),
  fHistClustersEnergy(0),
  fHistEPcorrelation(0),
  fHistJetsEnergy(0),
  fHistJetsNE(0),
  fHistJetsNEF(0),
  fHistJetsZ(0),
  fHistTrPhiEta(0),
  fHistClusPhiEta(0),
  fHistJetPhiEta(0),
  fHistMaxTrgCluster(0),
  Ptbins(100),
  Ptlow(0),
  Ptup(50),
  Ebins(100),
  Elow(0),
  Eup(50)
{
  // Standard constructor.

  if (!name)
    return;

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

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

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

//________________________________________________________________________
void AliEmcalIsolatedPhotonsTask::UserCreateOutputObjects()
{
   // Create histograms
  // Called once (on the worker node)
  
  fOutput = new TList();
  fOutput->SetOwner();  // IMPORTANT!
  
  AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
  
  if (handler && handler->InheritsFrom("AliESDInputHandler")) {

    if (!fESDTrackCuts) fESDTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
    
    if (!fESDTrackCuts) {
      AliFatal("Invalid pointer to ESDtrackCuts");
    }
    
    if (fESDTrackCuts) {
      fESDTrackCuts->DefineHistograms(kRed);
      fOutput->Add(fESDTrackCuts);
    }
  }
    
  fHistTracksPt = new TH1F("fHistTracksPt","P_{T} spectrum of reconstructed tracks", Ptbins, Ptlow, Ptup);
  fHistTracksPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
  fHistTracksPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistTracksPt);
	
  fHistClustersEnergy = new TH1F("fHistClustersEnergy","Energy spectrum of clusters", Ebins, Elow, Eup);
  fHistClustersEnergy->GetXaxis()->SetTitle("E [GeV]");
  fHistClustersEnergy->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistClustersEnergy);

  fHistEPcorrelation = new TH2F("fHistEPcorrelation","Energy-momentum correlation", Ptbins, Ptlow, Ptup, Ebins, Elow, Eup);
  fHistEPcorrelation->GetXaxis()->SetTitle("P [GeV/c]");
  fHistEPcorrelation->GetYaxis()->SetTitle("E [GeV]");
  fOutput->Add(fHistEPcorrelation);
  
  fHistJetsEnergy = new TH1F("fHistJetsEnergy","Energy spectrum of jets", Ebins, Elow, Eup);
  fHistJetsEnergy->GetXaxis()->SetTitle("E [GeV]");
  fHistJetsEnergy->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsEnergy);
  
  fHistJetsNE = new TH1F("fHistJetsNE","Neutral energy spectrum of jets", Ebins, Elow, Eup);
  fHistJetsNE->GetXaxis()->SetTitle("E [GeV]");
  fHistJetsNE->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsNE);
  
  fHistJetsNEF = new TH1F("fHistJetsNEF","Jets neutral energy fraction", Ebins, 0, 1.2);
  fHistJetsNEF->GetXaxis()->SetTitle("E [GeV]");
  fHistJetsNEF->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsNEF);

  fHistJetsZ = new TH1F("fHistJetsZ","Z of jet constituents", Ebins, 0, 1.2);
  fHistJetsZ->GetXaxis()->SetTitle("Z");
  fHistJetsZ->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsZ);

  fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 20, -2, 2, 32, 0, 6.4);
  fHistTrPhiEta->GetXaxis()->SetTitle("Eta");
  fHistTrPhiEta->GetYaxis()->SetTitle("Phi");
  fOutput->Add(fHistTrPhiEta);

  fHistClusPhiEta = new TH2F("fHistClusPhiEta","Phi-Eta distribution of clusters", 20, -2, 2, 32, 0, 6.4);
  fHistClusPhiEta->GetXaxis()->SetTitle("Eta");
  fHistClusPhiEta->GetYaxis()->SetTitle("Phi");
  fOutput->Add(fHistClusPhiEta);

  fHistJetPhiEta = new TH2F("fHistJetPhiEta","Phi-Eta distribution of jets", 20, -2, 2, 32, 0, 6.4);
  fHistJetPhiEta->GetXaxis()->SetTitle("Eta");
  fHistJetPhiEta->GetYaxis()->SetTitle("Phi");
  fOutput->Add(fHistJetPhiEta);

  fHistMaxTrgCluster = new TH1F("fHistMaxTrgCluster","Energy distribution of max trigger clusters", Ebins, Elow, Eup);
  fHistMaxTrgCluster->GetXaxis()->SetTitle("E [GeV]");
  fHistMaxTrgCluster->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistMaxTrgCluster);

  for (Int_t i = 0; i < 4; i++) {
    TString histnamephi("fHistTrackPhi_");
    histnamephi += i;
    fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 32, 0, 6.4);
    fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
    fOutput->Add(fHistTrackPhi[i]);

    TString histnameeta("fHistTrackEta_");
    histnameeta += i;
    fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 40, -2, 2);
    fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
    fOutput->Add(fHistTrackEta[i]);
  }

  fHistTrackPhi[0]->SetLineColor(kRed);
  fHistTrackEta[0]->SetLineColor(kRed);
  fHistTrackPhi[1]->SetLineColor(kBlue);
  fHistTrackEta[1]->SetLineColor(kBlue);
  fHistTrackPhi[2]->SetLineColor(kGreen);
  fHistTrackEta[2]->SetLineColor(kGreen);
  fHistTrackPhi[3]->SetLineColor(kBlack);
  fHistTrackEta[3]->SetLineColor(kBlack);

  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));
  fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
  fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
  fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));

  if (!fTracks) {
    AliWarning(Form("Could not retrieve tracks!")); 
  }
  if (!fCaloClusters) {
    AliWarning(Form("Could not retrieve clusters!")); 
  }
  if (!fJets) {
    AliWarning(Form("Could not retrieve jets!")); 
  }
  if (!fTrgClusters) {
    AliWarning(Form("Could not retrieve trigger clusters!")); 
  }
}

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()
{
  
  // Cluster loop
  Int_t nclusters =  GetNumberOfCaloClusters();
  //cout << nclusters << " clusters" << endl;
  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;

    fHistClustersEnergy->Fill(cluster->E());

    Float_t pos[3];
    cluster->GetPosition(pos);
    TVector3 clusVec(pos);
    fHistClusPhiEta->Fill(clusVec.Eta(), clusVec.Phi());
  } //cluster loop 
  

  // Track loop 
  Int_t ntracks = GetNumberOfTracks();
  //cout << ntracks << " tracks" << endl;
  for(Int_t i = 0; i < ntracks; i++) {
    //cout << "track n. " << i << endl;
    AliVTrack* track = GetTrack(i); // pointer to reconstructed to track          
    if(!track) {
      AliError(Form("ERROR: Could not retrieve esdtrack %d",i)); 
      continue; 
    }
    
    if (!AcceptTrack(track)) continue;

    fHistTracksPt->Fill(track->Pt());

    Int_t clId = track->GetEMCALcluster();
    if (clId > -1) {
      AliVCluster* cluster = GetCaloCluster(clId);
      if (cluster)
	fHistEPcorrelation->Fill(track->P(),cluster->E());
    } 

    Float_t eta,phi;
    Int_t label;

    if(track->InheritsFrom("AliESDtrack")) {
      AliESDtrack *esdtrack = dynamic_cast<AliESDtrack*>(track);
      eta = esdtrack->Eta();
      phi = esdtrack->Phi();
      label = esdtrack->GetLabel();
    }
    else if (track->InheritsFrom("AliAODTrack")) {
      AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
      eta = aodtrack->Eta();
      phi = aodtrack->Phi();
      label = aodtrack->GetLabel();
    }
    else if (track->InheritsFrom("AliPicoTrack")) {
      AliPicoTrack *picotrack = dynamic_cast<AliPicoTrack*>(track);
      eta = picotrack->Eta();
      phi = picotrack->Phi();
      label = picotrack->GetLabel();
    }
    else {
      AliWarning("Track type not recognized! Will not plot phi\eta distributions!");
      continue;
    }

    fHistTrPhiEta->Fill(eta, phi);
    fHistTrackEta[3]->Fill(eta);
    fHistTrackPhi[3]->Fill(phi);

    if (label >= 0 && label < 3) {
      fHistTrackEta[label]->Fill(eta);
      fHistTrackPhi[label]->Fill(phi);
    }
    else {
      AliWarning(Form("Track label %d not recognized!",label));
    }
  }

  // Jet loop
  Int_t njets =  GetNumberOfJets();
  //cout << njets << " jets" << endl;
  for (Int_t ij = 0; ij < njets; ij++) {
    AliEmcalJet* jet = GetJet(ij);
    if (!jet) {
      printf("ERROR: Could not receive jet %d\n", ij);
      continue;
    }  
    
    fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
    fHistJetsEnergy->Fill(jet->E());
    fHistJetsNEF->Fill(jet->NEF());
    fHistJetsNE->Fill(jet->E() * jet->NEF());
    //if (jet->E() <= 0)
    //continue;
    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
      Int_t trackid = jet->TrackAt(it);
      AliVTrack *track = GetTrack(trackid);
      if (track)
	fHistJetsZ->Fill(track->Pt() / jet->E());
    }
    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
      Int_t clusterid = jet->ClusterAt(ic);
      AliVCluster *cluster = GetCaloCluster(clusterid);
      if (cluster)
	fHistJetsZ->Fill(cluster->E() / jet->E());
    }
  } //jet loop 

  Int_t ntrgclusters =  GetNumberOfTrgClusters();
  Float_t maxe = 0;
  //cout << ntrgclusters << " clusters" << endl;
  for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
    AliVCluster* cluster = GetTrgCluster(iClusters);
    if (!cluster) {
      printf("ERROR: Could not receive cluster %d\n", iClusters);
      continue;
    }  
    
    if (!(cluster->IsEMCAL())) continue;

    if (cluster->E() > maxe)
      maxe = cluster->E();

  } //cluster loop 
  fHistMaxTrgCluster->Fill(maxe);
  
}

//________________________________________________________________________
Bool_t AliEmcalIsolatedPhotonsTask::AcceptTrack(AliVTrack *track)
{
  if (!strcmp(track->ClassName(), "AliESDTrack") && fESDTrackCuts) {
    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack*>(track);
    if(esdtrack) {
      return fESDTrackCuts->AcceptTrack(esdtrack);
    }
    //else {
    //cout << "no esdtrack!" << endl;
    //}
  }
  else if (!strcmp(track->ClassName(), "AliAODTrack")) {
    AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
    if (aodtrack) {
      //cout << "filter bit = " << fFilterBit << ", filter map = " << aodtrack->GetFilterMap() << endl;
      return aodtrack->TestFilterBit(fFilterBit);
      
    }
  }
  return 1;
}

//________________________________________________________________________
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"
	13	#include "AliESDCaloCluster.h"
	14	#include "AliESDtrack.h"
	15	#include "AliEmcalJet.h"
	16	#include "AliFJWrapper.h"
	17	#include "AliAODTrack.h"
e82e282c	18	#include "AliESDtrack.h"
25283b37	19	#include "AliESDtrackCuts.h"
	20	#include "AliEmcalJet.h"
	21	#include "AliVEventHandler.h"
35789a2d	22	#include "AliPicoTrack.h"
25283b37	23
	24	#include "AliEmcalIsolatedPhotonsTask.h"
	25
	26	ClassImp(AliEmcalIsolatedPhotonsTask)
	27
	28	//________________________________________________________________________
	29	AliEmcalIsolatedPhotonsTask::AliEmcalIsolatedPhotonsTask() :
	30	AliAnalysisTaskSE("AliEmcalIsolatedPhotonsTask"),
	31	fOutput(0),
	32	fTracksName("Tracks"),
	33	fCaloName("CaloClusters"),
	34	fJetsName("Jets"),
e82e282c	35	fTrgClusName("ClustersL1GAMMAFEE"),
25283b37	36	fESDTrackCuts(0),
	37	fFilterBit(16),
	38	fTracks(0),
	39	fCaloClusters(0),
	40	fJets(0),
e82e282c	41	fTrgClusters(0),
25283b37	42	fHistTracksPt(0),
	43	fHistClustersEnergy(0),
	44	fHistEPcorrelation(0),
	45	fHistJetsEnergy(0),
	46	fHistJetsNE(0),
	47	fHistJetsNEF(0),
	48	fHistJetsZ(0),
	49	fHistTrPhiEta(0),
	50	fHistClusPhiEta(0),
	51	fHistJetPhiEta(0),
e82e282c	52	fHistMaxTrgCluster(0),
25283b37	53	Ptbins(100),
	54	Ptlow(0),
	55	Ptup(50),
	56	Ebins(100),
	57	Elow(0),
	58	Eup(50)
	59	{
	60	// Default constructor.
	61	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
	62	}
	63
	64	//________________________________________________________________________
	65	AliEmcalIsolatedPhotonsTask::AliEmcalIsolatedPhotonsTask(const char *name) :
	66	AliAnalysisTaskSE("AliEmcalIsolatedPhotonsTask"),
	67	fOutput(0),
	68	fTracksName("Tracks"),
	69	fCaloName("CaloClusters"),
	70	fJetsName("Jets"),
e82e282c	71	fTrgClusName("ClustersL1GAMMAFEE"),
25283b37	72	fESDTrackCuts(0),
	73	fFilterBit(16),
	74	fTracks(0),
	75	fCaloClusters(0),
	76	fJets(0),
e82e282c	77	fTrgClusters(0),
25283b37	78	fHistTracksPt(0),
	79	fHistClustersEnergy(0),
	80	fHistEPcorrelation(0),
	81	fHistJetsEnergy(0),
	82	fHistJetsNE(0),
	83	fHistJetsNEF(0),
	84	fHistJetsZ(0),
	85	fHistTrPhiEta(0),
	86	fHistClusPhiEta(0),
	87	fHistJetPhiEta(0),
e82e282c	88	fHistMaxTrgCluster(0),
25283b37	89	Ptbins(100),
	90	Ptlow(0),
	91	Ptup(50),
	92	Ebins(100),
	93	Elow(0),
	94	Eup(50)
	95	{
	96	// Standard constructor.
	97
	98	if (!name)
	99	return;
	100
	101	SetName(name);
	102	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
	103
	104	DefineInput(0,TChain::Class());
	105	DefineOutput(1,TList::Class());
	106	}
	107
	108	//________________________________________________________________________
	109	AliEmcalIsolatedPhotonsTask::~AliEmcalIsolatedPhotonsTask()
	110	{
	111	// Destructor
	112	}
	113
	114	//________________________________________________________________________
	115	void AliEmcalIsolatedPhotonsTask::UserCreateOutputObjects()
	116	{
	117	// Create histograms
	118	// Called once (on the worker node)
	119
	120	fOutput = new TList();
	121	fOutput->SetOwner(); // IMPORTANT!
	122
	123	AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
	124
e82e282c	125	if (handler && handler->InheritsFrom("AliESDInputHandler")) {
25283b37	126
	127	if (!fESDTrackCuts) fESDTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
	128
	129	if (!fESDTrackCuts) {
	130	AliFatal("Invalid pointer to ESDtrackCuts");
	131	}
	132
	133	if (fESDTrackCuts) {
	134	fESDTrackCuts->DefineHistograms(kRed);
	135	fOutput->Add(fESDTrackCuts);
	136	}
	137	}
25283b37	138
	139	fHistTracksPt = new TH1F("fHistTracksPt","P_{T} spectrum of reconstructed tracks", Ptbins, Ptlow, Ptup);
	140	fHistTracksPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
	141	fHistTracksPt->GetYaxis()->SetTitle("counts");
	142	fOutput->Add(fHistTracksPt);
	143
	144	fHistClustersEnergy = new TH1F("fHistClustersEnergy","Energy spectrum of clusters", Ebins, Elow, Eup);
	145	fHistClustersEnergy->GetXaxis()->SetTitle("E [GeV]");
	146	fHistClustersEnergy->GetYaxis()->SetTitle("counts");
	147	fOutput->Add(fHistClustersEnergy);
	148
	149	fHistEPcorrelation = new TH2F("fHistEPcorrelation","Energy-momentum correlation", Ptbins, Ptlow, Ptup, Ebins, Elow, Eup);
	150	fHistEPcorrelation->GetXaxis()->SetTitle("P [GeV/c]");
	151	fHistEPcorrelation->GetYaxis()->SetTitle("E [GeV]");
	152	fOutput->Add(fHistEPcorrelation);
	153
	154	fHistJetsEnergy = new TH1F("fHistJetsEnergy","Energy spectrum of jets", Ebins, Elow, Eup);
	155	fHistJetsEnergy->GetXaxis()->SetTitle("E [GeV]");
	156	fHistJetsEnergy->GetYaxis()->SetTitle("counts");
	157	fOutput->Add(fHistJetsEnergy);
	158
	159	fHistJetsNE = new TH1F("fHistJetsNE","Neutral energy spectrum of jets", Ebins, Elow, Eup);
	160	fHistJetsNE->GetXaxis()->SetTitle("E [GeV]");
	161	fHistJetsNE->GetYaxis()->SetTitle("counts");
	162	fOutput->Add(fHistJetsNE);
	163
	164	fHistJetsNEF = new TH1F("fHistJetsNEF","Jets neutral energy fraction", Ebins, 0, 1.2);
	165	fHistJetsNEF->GetXaxis()->SetTitle("E [GeV]");
	166	fHistJetsNEF->GetYaxis()->SetTitle("counts");
	167	fOutput->Add(fHistJetsNEF);
	168
	169	fHistJetsZ = new TH1F("fHistJetsZ","Z of jet constituents", Ebins, 0, 1.2);
	170	fHistJetsZ->GetXaxis()->SetTitle("Z");
	171	fHistJetsZ->GetYaxis()->SetTitle("counts");
	172	fOutput->Add(fHistJetsZ);
	173
	174	fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 20, -2, 2, 32, 0, 6.4);
	175	fHistTrPhiEta->GetXaxis()->SetTitle("Eta");
	176	fHistTrPhiEta->GetYaxis()->SetTitle("Phi");
	177	fOutput->Add(fHistTrPhiEta);
	178
	179	fHistClusPhiEta = new TH2F("fHistClusPhiEta","Phi-Eta distribution of clusters", 20, -2, 2, 32, 0, 6.4);
	180	fHistClusPhiEta->GetXaxis()->SetTitle("Eta");
	181	fHistClusPhiEta->GetYaxis()->SetTitle("Phi");
	182	fOutput->Add(fHistClusPhiEta);
	183
	184	fHistJetPhiEta = new TH2F("fHistJetPhiEta","Phi-Eta distribution of jets", 20, -2, 2, 32, 0, 6.4);
	185	fHistJetPhiEta->GetXaxis()->SetTitle("Eta");
	186	fHistJetPhiEta->GetYaxis()->SetTitle("Phi");
	187	fOutput->Add(fHistJetPhiEta);
e82e282c	188
	189	fHistMaxTrgCluster = new TH1F("fHistMaxTrgCluster","Energy distribution of max trigger clusters", Ebins, Elow, Eup);
	190	fHistMaxTrgCluster->GetXaxis()->SetTitle("E [GeV]");
	191	fHistMaxTrgCluster->GetYaxis()->SetTitle("counts");
	192	fOutput->Add(fHistMaxTrgCluster);
	193
680855e5	194	for (Int_t i = 0; i < 4; i++) {
e82e282c	195	TString histnamephi("fHistTrackPhi_");
	196	histnamephi += i;
	197	fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 32, 0, 6.4);
	198	fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
	199	fOutput->Add(fHistTrackPhi[i]);
	200
	201	TString histnameeta("fHistTrackEta_");
	202	histnameeta += i;
	203	fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 40, -2, 2);
	204	fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
	205	fOutput->Add(fHistTrackEta[i]);
	206	}
	207
	208	fHistTrackPhi[0]->SetLineColor(kRed);
	209	fHistTrackEta[0]->SetLineColor(kRed);
	210	fHistTrackPhi[1]->SetLineColor(kBlue);
	211	fHistTrackEta[1]->SetLineColor(kBlue);
	212	fHistTrackPhi[2]->SetLineColor(kGreen);
	213	fHistTrackEta[2]->SetLineColor(kGreen);
680855e5	214	fHistTrackPhi[3]->SetLineColor(kBlack);
680855e5	215	fHistTrackEta[3]->SetLineColor(kBlack);
e82e282c	216
25283b37	217	PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
	218	}
	219
	220	void AliEmcalIsolatedPhotonsTask::RetrieveEventObjects()
	221	{
	222	fCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
	223	fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
	224	fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
e82e282c	225	fTrgClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
25283b37	226
25283b37	227	if (!fTracks) {
e82e282c	228	AliWarning(Form("Could not retrieve tracks!"));
25283b37	229	}
25283b37	230	if (!fCaloClusters) {
e82e282c	231	AliWarning(Form("Could not retrieve clusters!"));
25283b37	232	}
25283b37	233	if (!fJets) {
e82e282c	234	AliWarning(Form("Could not retrieve jets!"));
	235	}
	236	if (!fTrgClusters) {
	237	AliWarning(Form("Could not retrieve trigger clusters!"));
25283b37	238	}
	239	}
	240
	241	AliVTrack* AliEmcalIsolatedPhotonsTask::GetTrack(const Int_t i) const
	242	{
e82e282c	243	if (fTracks)
	244	return dynamic_cast<AliVTrack*>(fTracks->At(i));
	245	else
	246	return 0;
25283b37	247	}
	248
	249	Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfTracks() const
	250	{
e82e282c	251	if (fTracks)
	252	return fTracks->GetEntriesFast();
	253	else
	254	return 0;
25283b37	255	}
	256
	257	AliVCluster* AliEmcalIsolatedPhotonsTask::GetCaloCluster(const Int_t i) const
e82e282c	258	{
	259	if (fCaloClusters)
	260	return dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
	261	else
	262	return 0;
25283b37	263	}
	264
	265	Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfCaloClusters() const
e82e282c	266	{
	267	if (fCaloClusters)
	268	return fCaloClusters->GetEntriesFast();
	269	else
	270	return 0;
25283b37	271	}
	272
	273	AliEmcalJet* AliEmcalIsolatedPhotonsTask::GetJet(const Int_t i) const
	274	{
e82e282c	275	if (fJets)
	276	return dynamic_cast<AliEmcalJet*>(fJets->At(i));
	277	else
	278	return 0;
25283b37	279	}
	280
	281	Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfJets() const
	282	{
e82e282c	283	if (fJets)
	284	return fJets->GetEntriesFast();
	285	else
	286	return 0;
	287	}
	288
	289	AliVCluster* AliEmcalIsolatedPhotonsTask::GetTrgCluster(const Int_t i) const
	290	{
	291	if (fTrgClusters)
	292	return dynamic_cast<AliVCluster*>(fTrgClusters->At(i));
	293	else
	294	return 0;
	295	}
	296
	297	Int_t AliEmcalIsolatedPhotonsTask::GetNumberOfTrgClusters() const
	298	{
	299	if (fTrgClusters)
	300	return fTrgClusters->GetEntriesFast();
	301	else
	302	return 0;
25283b37	303	}
	304
	305	void AliEmcalIsolatedPhotonsTask::FillHistograms()
	306	{
	307
	308	// Cluster loop
	309	Int_t nclusters = GetNumberOfCaloClusters();
	310	//cout << nclusters << " clusters" << endl;
	311	for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
	312	AliVCluster* cluster = GetCaloCluster(iClusters);
	313	if (!cluster) {
	314	printf("ERROR: Could not receive cluster %d\n", iClusters);
	315	continue;
	316	}
	317
	318	if (!(cluster->IsEMCAL())) continue;
	319
	320	fHistClustersEnergy->Fill(cluster->E());
	321
	322	Float_t pos[3];
	323	cluster->GetPosition(pos);
	324	TVector3 clusVec(pos);
	325	fHistClusPhiEta->Fill(clusVec.Eta(), clusVec.Phi());
	326	} //cluster loop
	327
	328
	329	// Track loop
	330	Int_t ntracks = GetNumberOfTracks();
	331	//cout << ntracks << " tracks" << endl;
	332	for(Int_t i = 0; i < ntracks; i++) {
	333	//cout << "track n. " << i << endl;
	334	AliVTrack* track = GetTrack(i); // pointer to reconstructed to track
	335	if(!track) {
	336	AliError(Form("ERROR: Could not retrieve esdtrack %d",i));
	337	continue;
	338	}
	339
	340	if (!AcceptTrack(track)) continue;
	341
	342	fHistTracksPt->Fill(track->Pt());
e82e282c	343
	344	Int_t clId = track->GetEMCALcluster();
	345	if (clId > -1) {
	346	AliVCluster* cluster = GetCaloCluster(clId);
	347	if (cluster)
	348	fHistEPcorrelation->Fill(track->P(),cluster->E());
	349	}
	350
	351	Float_t eta,phi;
	352	Int_t label;
	353
	354	if(track->InheritsFrom("AliESDtrack")) {
	355	AliESDtrack esdtrack = dynamic_cast<AliESDtrack>(track);
	356	eta = esdtrack->Eta();
	357	phi = esdtrack->Phi();
	358	label = esdtrack->GetLabel();
35789a2d	359	}
	360	else if (track->InheritsFrom("AliAODTrack")) {
	361	AliAODTrack aodtrack = dynamic_cast<AliAODTrack>(track);
e82e282c	362	eta = aodtrack->Eta();
	363	phi = aodtrack->Phi();
	364	label = aodtrack->GetLabel();
35789a2d	365	}
	366	else if (track->InheritsFrom("AliPicoTrack")) {
	367	AliPicoTrack picotrack = dynamic_cast<AliPicoTrack>(track);
e82e282c	368	eta = picotrack->Eta();
	369	phi = picotrack->Phi();
	370	label = picotrack->GetLabel();
	371	}
	372	else {
	373	AliWarning("Track type not recognized! Will not plot phi\eta distributions!");
	374	continue;
35789a2d	375	}
25283b37	376
e82e282c	377	fHistTrPhiEta->Fill(eta, phi);
680855e5	378	fHistTrackEta[3]->Fill(eta);
680855e5	379	fHistTrackPhi[3]->Fill(phi);
e82e282c	380
	381	if (label >= 0 && label < 3) {
	382	fHistTrackEta[label]->Fill(eta);
	383	fHistTrackPhi[label]->Fill(phi);
	384	}
	385	else {
	386	AliWarning(Form("Track label %d not recognized!",label));
	387	}
25283b37	388	}
	389
	390	// Jet loop
	391	Int_t njets = GetNumberOfJets();
	392	//cout << njets << " jets" << endl;
	393	for (Int_t ij = 0; ij < njets; ij++) {
	394	AliEmcalJet* jet = GetJet(ij);
	395	if (!jet) {
	396	printf("ERROR: Could not receive jet %d\n", ij);
	397	continue;
	398	}
	399
	400	fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
	401	fHistJetsEnergy->Fill(jet->E());
	402	fHistJetsNEF->Fill(jet->NEF());
	403	fHistJetsNE->Fill(jet->E() * jet->NEF());
35789a2d	404	//if (jet->E() <= 0)
	405	//continue;
	406	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
	407	Int_t trackid = jet->TrackAt(it);
	408	AliVTrack *track = GetTrack(trackid);
	409	if (track)
	410	fHistJetsZ->Fill(track->Pt() / jet->E());
	411	}
	412	for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
	413	Int_t clusterid = jet->ClusterAt(ic);
	414	AliVCluster *cluster = GetCaloCluster(clusterid);
	415	if (cluster)
	416	fHistJetsZ->Fill(cluster->E() / jet->E());
	417	}
25283b37	418	} //jet loop
e82e282c	419
	420	Int_t ntrgclusters = GetNumberOfTrgClusters();
	421	Float_t maxe = 0;
	422	//cout << ntrgclusters << " clusters" << endl;
	423	for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
	424	AliVCluster* cluster = GetTrgCluster(iClusters);
	425	if (!cluster) {
	426	printf("ERROR: Could not receive cluster %d\n", iClusters);
	427	continue;
	428	}
	429
	430	if (!(cluster->IsEMCAL())) continue;
	431
	432	if (cluster->E() > maxe)
	433	maxe = cluster->E();
	434
	435	} //cluster loop
	436	fHistMaxTrgCluster->Fill(maxe);
	437
25283b37	438	}
	439
	440	//________________________________________________________________________
	441	Bool_t AliEmcalIsolatedPhotonsTask::AcceptTrack(AliVTrack *track)
	442	{
25283b37	443	if (!strcmp(track->ClassName(), "AliESDTrack") && fESDTrackCuts) {
	444	AliESDtrack esdtrack = dynamic_cast<AliESDtrack>(track);
	445	if(esdtrack) {
	446	return fESDTrackCuts->AcceptTrack(esdtrack);
	447	}
	448	//else {
	449	//cout << "no esdtrack!" << endl;
	450	//}
	451	}
	452	else if (!strcmp(track->ClassName(), "AliAODTrack")) {
	453	AliAODTrack aodtrack = dynamic_cast<AliAODTrack>(track);
	454	if (aodtrack) {
	455	//cout << "filter bit = " << fFilterBit << ", filter map = " << aodtrack->GetFilterMap() << endl;
	456	return aodtrack->TestFilterBit(fFilterBit);
	457
	458	}
	459	}
	460	return 1;
	461	}
	462
	463	//________________________________________________________________________
	464	void AliEmcalIsolatedPhotonsTask::UserExec(Option_t *)
	465	{
	466	// Main loop, called for each event.
	467	// Add jets to event if not yet there
	468
	469	RetrieveEventObjects();
	470
	471	FillHistograms();
	472
	473	// information for this iteration of the UserExec in the container
	474	PostData(1, fOutput);
	475	}
	476
	477	//________________________________________________________________________
	478	void AliEmcalIsolatedPhotonsTask::Terminate(Option_t *)
	479	{
	480	// Called once at the end of the analysis.
	481	}