[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskEmcalJetSample.cxx

// $Id$
//
// Jet sample analysis task.
//
// Author: S.Aiola

#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TList.h>
#include <TLorentzVector.h>

#include "AliVCluster.h"
#include "AliAODCaloCluster.h"
#include "AliESDCaloCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliJetContainer.h"
#include "AliParticleContainer.h"
#include "AliClusterContainer.h"
#include "AliPicoTrack.h"

#include "AliAnalysisTaskEmcalJetSample.h"

ClassImp(AliAnalysisTaskEmcalJetSample)

//________________________________________________________________________
AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE),
  fHistTracksPt(0),
  fHistClustersPt(0),
  fHistLeadingJetPt(0),
  fHistJetsPhiEta(0),
  fHistJetsPtArea(0),
  fHistJetsPtLeadHad(0),
  fHistJetsCorrPtArea(0),
  fHistPtDEtaDPhiTrackClus(0),
  fHistPtDEtaDPhiClusTrack(0),
  fHistClustDx(0),
  fHistClustDz(0),
  fJetsCont(0),
  fTracksCont(0),
  fCaloClustersCont(0)

{
  // Default constructor.

  fHistTracksPt       = new TH1*[fNcentBins];
  fHistClustersPt     = new TH1*[fNcentBins];
  fHistLeadingJetPt   = new TH1*[fNcentBins];
  fHistJetsPhiEta     = new TH2*[fNcentBins];
  fHistJetsPtArea     = new TH2*[fNcentBins];
  fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  fHistJetsCorrPtArea = new TH2*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fHistTracksPt[i] = 0;
    fHistClustersPt[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHistJetsPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsPtLeadHad[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),
  fHistTracksPt(0),
  fHistClustersPt(0),
  fHistLeadingJetPt(0),
  fHistJetsPhiEta(0),
  fHistJetsPtArea(0),
  fHistJetsPtLeadHad(0),
  fHistJetsCorrPtArea(0),
  fHistPtDEtaDPhiTrackClus(0),
  fHistPtDEtaDPhiClusTrack(0),
  fHistClustDx(0),
  fHistClustDz(0),
  fJetsCont(0),
  fTracksCont(0),
  fCaloClustersCont(0)
{
  // Standard constructor.

  fHistTracksPt       = new TH1*[fNcentBins];
  fHistClustersPt     = new TH1*[fNcentBins];
  fHistLeadingJetPt   = new TH1*[fNcentBins];
  fHistJetsPhiEta     = new TH2*[fNcentBins];
  fHistJetsPtArea     = new TH2*[fNcentBins];
  fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  fHistJetsCorrPtArea = new TH2*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fHistTracksPt[i] = 0;
    fHistClustersPt[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHistJetsPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsPtLeadHad[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetSample::~AliAnalysisTaskEmcalJetSample()
{
  // Destructor.
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetSample::UserCreateOutputObjects()
{
  // Create user output.

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  fJetsCont           = GetJetContainer(0);
  if(fJetsCont) { //get particles and clusters connected to jets
    fTracksCont       = fJetsCont->GetParticleContainer();
    fCaloClustersCont = fJetsCont->GetClusterContainer();
  } else {        //no jets, just analysis tracks and clusters
    fTracksCont       = GetParticleContainer(0);
    fCaloClustersCont = GetClusterContainer(0);
  }
  if(fTracksCont) fTracksCont->SetClassName("AliVTrack");
  if(fCaloClustersCont) fCaloClustersCont->SetClassName("AliVCluster");

  TString histname;

  for (Int_t i = 0; i < fNcentBins; i++) {
    if (fParticleCollArray.GetEntriesFast()>0) {
      histname = "fHistTracksPt_";
      histname += i;
      fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
      fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
      fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistTracksPt[i]);
    }

    if (fClusterCollArray.GetEntriesFast()>0) {
      histname = "fHistClustersPt_";
      histname += i;
      fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
      fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
      fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistClustersPt[i]);
    }

    if (fJetCollArray.GetEntriesFast()>0) {
      histname = "fHistLeadingJetPt_";
      histname += i;
      fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
      fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
      fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistLeadingJetPt[i]);
      
      histname = "fHistJetsPhiEta_";
      histname += i;
      fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
      fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
      fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
      fOutput->Add(fHistJetsPhiEta[i]);
      
      histname = "fHistJetsPtArea_";
      histname += i;
      fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
      fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
      fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
      fOutput->Add(fHistJetsPtArea[i]);

      histname = "fHistJetsPtLeadHad_";
      histname += i;
      fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
      fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
      fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
      fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
      fOutput->Add(fHistJetsPtLeadHad[i]);
    
      if (!(GetJetContainer()->GetRhoName().IsNull())) {
	histname = "fHistJetsCorrPtArea_";
	histname += i;
	fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
	fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
	fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
	fOutput->Add(fHistJetsCorrPtArea[i]);
      }
    }
  }

  histname = "fHistPtDEtaDPhiTrackClus";
  fHistPtDEtaDPhiTrackClus = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{track};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
  fOutput->Add(fHistPtDEtaDPhiTrackClus);

  histname = "fHistPtDEtaDPhiClusTrack";
  fHistPtDEtaDPhiClusTrack = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{clus};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
  fOutput->Add(fHistPtDEtaDPhiClusTrack);

  fHistClustDx = new TH1F("fHistClustDx","fHistClustDx;Dx",1000,0.,1.);
  fOutput->Add(fHistClustDx);

  fHistClustDz = new TH1F("fHistClustDz","fHistClustDz;Dz",1000,0.,1.);
  fOutput->Add(fHistClustDz);

  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetSample::FillHistograms()
{
  // Fill histograms.

  if (fTracksCont) {
    AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0)); 
    while(track) {
      fHistTracksPt[fCentBin]->Fill(track->Pt()); 
      track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
    }
  }

  if (fCaloClustersCont) {
    AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0); 
    while(cluster) {
      TLorentzVector nPart;
      cluster->GetMomentum(nPart, fVertex);
      fHistClustersPt[fCentBin]->Fill(nPart.Pt());
      Double_t dx = cluster->GetTrackDx();
      Double_t dz = cluster->GetTrackDz();
      fHistClustDx->Fill(dx);
      fHistClustDz->Fill(dz);
      cluster = fCaloClustersCont->GetNextAcceptCluster();
    }
  }

  if (fJetsCont) {
    AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0); 
    while(jet) {

      fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
      fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());

      Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
      fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);

      if (fHistJetsCorrPtArea[fCentBin]) {
	Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
	fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
      }
      jet = fJetsCont->GetNextAcceptJet(); 
    }
    
    jet = fJetsCont->GetLeadingJet();
    if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
  }

  CheckClusTrackMatching();

  return kTRUE;
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetSample::CheckClusTrackMatching()
{
  
  if(!fTracksCont || !fCaloClustersCont)
    return;

  Double_t deta = 999;
  Double_t dphi = 999;

  //Get closest cluster to track
  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0)); 
  while(track) {
    //Get matched cluster
    Int_t emc1 = track->GetEMCALcluster();
    if(fCaloClustersCont && emc1>=0) {
      AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
      if(clusMatch) {
	AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
	fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
      }
    }
    track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
  }
  
  //Get closest track to cluster
  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0); 
  while(cluster) {
    TLorentzVector nPart;
    cluster->GetMomentum(nPart, fVertex);
    fHistClustersPt[fCentBin]->Fill(nPart.Pt());
    
    //Get matched track
    AliVTrack *mt = NULL;      
    AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
    if(acl) {
      if(acl->GetNTracksMatched()>1)
	mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
    }
    else {
      AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
      Int_t im = ecl->GetTrackMatchedIndex();
      if(fTracksCont && im>=0) {
	mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
      }
    }
    if(mt) {
      AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
      fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
      
      /* //debugging
	 if(mt->IsEMCAL()) {
	 Int_t emc1 = mt->GetEMCALcluster();
	 Printf("current id: %d  emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
	 AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
	 AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
	 Printf("deta: %f dphi: %f",deta,dphi);
	 }
      */
    }
    cluster = fCaloClustersCont->GetNextAcceptCluster();
  }
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetSample::ExecOnce() {

  AliAnalysisTaskEmcalJet::ExecOnce();

  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;

}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetSample::Run()
{
  // Run analysis code here, if needed. It will be executed before FillHistograms().

  return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetSample::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
8628b70c	1	// $Id$
9f52c61f	2	//
	3	// Jet sample analysis task.
	4	//
	5	// Author: S.Aiola
	6
	7	#include <TClonesArray.h>
	8	#include <TH1F.h>
	9	#include <TH2F.h>
0925db21	10	#include <TH3F.h>
9f52c61f	11	#include <TList.h>
	12	#include <TLorentzVector.h>
	13
	14	#include "AliVCluster.h"
4bb8c6df	15	#include "AliAODCaloCluster.h"
0925db21	16	#include "AliESDCaloCluster.h"
9f52c61f	17	#include "AliVTrack.h"
	18	#include "AliEmcalJet.h"
	19	#include "AliRhoParameter.h"
	20	#include "AliLog.h"
9239b066	21	#include "AliJetContainer.h"
8612dfc8	22	#include "AliParticleContainer.h"
8612dfc8	23	#include "AliClusterContainer.h"
0925db21	24	#include "AliPicoTrack.h"
9f52c61f	25
	26	#include "AliAnalysisTaskEmcalJetSample.h"
	27
	28	ClassImp(AliAnalysisTaskEmcalJetSample)
	29
	30	//________________________________________________________________________
	31	AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample() :
8612dfc8	32	AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE),
4bb8c6df	33	fHistTracksPt(0),
	34	fHistClustersPt(0),
	35	fHistLeadingJetPt(0),
	36	fHistJetsPhiEta(0),
	37	fHistJetsPtArea(0),
	38	fHistJetsPtLeadHad(0),
	39	fHistJetsCorrPtArea(0),
0925db21	40	fHistPtDEtaDPhiTrackClus(0),
0925db21	41	fHistPtDEtaDPhiClusTrack(0),
46b40b53	42	fHistClustDx(0),
46b40b53	43	fHistClustDz(0),
8612dfc8	44	fJetsCont(0),
	45	fTracksCont(0),
	46	fCaloClustersCont(0)
9f52c61f	47
	48	{
	49	// Default constructor.
	50
4bb8c6df	51	fHistTracksPt = new TH1*[fNcentBins];
	52	fHistClustersPt = new TH1*[fNcentBins];
	53	fHistLeadingJetPt = new TH1*[fNcentBins];
	54	fHistJetsPhiEta = new TH2*[fNcentBins];
	55	fHistJetsPtArea = new TH2*[fNcentBins];
	56	fHistJetsPtLeadHad = new TH2*[fNcentBins];
	57	fHistJetsCorrPtArea = new TH2*[fNcentBins];
	58
	59	for (Int_t i = 0; i < fNcentBins; i++) {
9f52c61f	60	fHistTracksPt[i] = 0;
	61	fHistClustersPt[i] = 0;
	62	fHistLeadingJetPt[i] = 0;
	63	fHistJetsPhiEta[i] = 0;
	64	fHistJetsPtArea[i] = 0;
	65	fHistJetsPtLeadHad[i] = 0;
	66	fHistJetsCorrPtArea[i] = 0;
	67	}
	68
	69	SetMakeGeneralHistograms(kTRUE);
	70	}
	71
	72	//________________________________________________________________________
	73	AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample(const char *name) :
8612dfc8	74	AliAnalysisTaskEmcalJet(name, kTRUE),
4bb8c6df	75	fHistTracksPt(0),
	76	fHistClustersPt(0),
	77	fHistLeadingJetPt(0),
	78	fHistJetsPhiEta(0),
	79	fHistJetsPtArea(0),
	80	fHistJetsPtLeadHad(0),
	81	fHistJetsCorrPtArea(0),
0925db21	82	fHistPtDEtaDPhiTrackClus(0),
0925db21	83	fHistPtDEtaDPhiClusTrack(0),
46b40b53	84	fHistClustDx(0),
46b40b53	85	fHistClustDz(0),
8612dfc8	86	fJetsCont(0),
	87	fTracksCont(0),
	88	fCaloClustersCont(0)
9f52c61f	89	{
	90	// Standard constructor.
	91
4bb8c6df	92	fHistTracksPt = new TH1*[fNcentBins];
	93	fHistClustersPt = new TH1*[fNcentBins];
	94	fHistLeadingJetPt = new TH1*[fNcentBins];
	95	fHistJetsPhiEta = new TH2*[fNcentBins];
	96	fHistJetsPtArea = new TH2*[fNcentBins];
	97	fHistJetsPtLeadHad = new TH2*[fNcentBins];
	98	fHistJetsCorrPtArea = new TH2*[fNcentBins];
9239b066	99
4bb8c6df	100	for (Int_t i = 0; i < fNcentBins; i++) {
9f52c61f	101	fHistTracksPt[i] = 0;
	102	fHistClustersPt[i] = 0;
	103	fHistLeadingJetPt[i] = 0;
	104	fHistJetsPhiEta[i] = 0;
	105	fHistJetsPtArea[i] = 0;
	106	fHistJetsPtLeadHad[i] = 0;
	107	fHistJetsCorrPtArea[i] = 0;
	108	}
	109
	110	SetMakeGeneralHistograms(kTRUE);
	111	}
	112
	113	//________________________________________________________________________
	114	AliAnalysisTaskEmcalJetSample::~AliAnalysisTaskEmcalJetSample()
	115	{
	116	// Destructor.
	117	}
	118
	119	//________________________________________________________________________
	120	void AliAnalysisTaskEmcalJetSample::UserCreateOutputObjects()
	121	{
	122	// Create user output.
	123
	124	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
	125
4bb8c6df	126	fJetsCont = GetJetContainer(0);
	127	if(fJetsCont) { //get particles and clusters connected to jets
	128	fTracksCont = fJetsCont->GetParticleContainer();
	129	fCaloClustersCont = fJetsCont->GetClusterContainer();
	130	} else { //no jets, just analysis tracks and clusters
	131	fTracksCont = GetParticleContainer(0);
	132	fCaloClustersCont = GetClusterContainer(0);
	133	}
3a0055e7	134	if(fTracksCont) fTracksCont->SetClassName("AliVTrack");
3a0055e7	135	if(fCaloClustersCont) fCaloClustersCont->SetClassName("AliVCluster");
8612dfc8	136
9f52c61f	137	TString histname;
9f52c61f	138
4bb8c6df	139	for (Int_t i = 0; i < fNcentBins; i++) {
9239b066	140	if (fParticleCollArray.GetEntriesFast()>0) {
9f52c61f	141	histname = "fHistTracksPt_";
	142	histname += i;
	143	fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
	144	fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
	145	fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
	146	fOutput->Add(fHistTracksPt[i]);
	147	}
	148
9239b066	149	if (fClusterCollArray.GetEntriesFast()>0) {
9f52c61f	150	histname = "fHistClustersPt_";
	151	histname += i;
	152	fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
	153	fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
	154	fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
	155	fOutput->Add(fHistClustersPt[i]);
	156	}
	157
9239b066	158	if (fJetCollArray.GetEntriesFast()>0) {
9f52c61f	159	histname = "fHistLeadingJetPt_";
	160	histname += i;
	161	fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
	162	fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
	163	fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
	164	fOutput->Add(fHistLeadingJetPt[i]);
	165
	166	histname = "fHistJetsPhiEta_";
	167	histname += i;
	168	fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
	169	fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
	170	fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
	171	fOutput->Add(fHistJetsPhiEta[i]);
	172
	173	histname = "fHistJetsPtArea_";
	174	histname += i;
	175	fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
	176	fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
	177	fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
	178	fOutput->Add(fHistJetsPtArea[i]);
	179
	180	histname = "fHistJetsPtLeadHad_";
	181	histname += i;
	182	fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
	183	fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
	184	fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
	185	fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
	186	fOutput->Add(fHistJetsPtLeadHad[i]);
	187
9239b066	188	if (!(GetJetContainer()->GetRhoName().IsNull())) {
9f52c61f	189	histname = "fHistJetsCorrPtArea_";
	190	histname += i;
	191	fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
	192	fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
	193	fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
	194	fOutput->Add(fHistJetsCorrPtArea[i]);
	195	}
	196	}
	197	}
0925db21	198
	199	histname = "fHistPtDEtaDPhiTrackClus";
	200	fHistPtDEtaDPhiTrackClus = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{track};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
	201	fOutput->Add(fHistPtDEtaDPhiTrackClus);
	202
	203	histname = "fHistPtDEtaDPhiClusTrack";
	204	fHistPtDEtaDPhiClusTrack = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{clus};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
	205	fOutput->Add(fHistPtDEtaDPhiClusTrack);
	206
46b40b53	207	fHistClustDx = new TH1F("fHistClustDx","fHistClustDx;Dx",1000,0.,1.);
	208	fOutput->Add(fHistClustDx);
	209
	210	fHistClustDz = new TH1F("fHistClustDz","fHistClustDz;Dz",1000,0.,1.);
	211	fOutput->Add(fHistClustDz);
	212
9f52c61f	213	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
	214	}
	215
	216	//________________________________________________________________________
	217	Bool_t AliAnalysisTaskEmcalJetSample::FillHistograms()
	218	{
	219	// Fill histograms.
	220
8612dfc8	221	if (fTracksCont) {
4bb8c6df	222	AliVTrack track = static_cast<AliVTrack>(fTracksCont->GetNextAcceptParticle(0));
8612dfc8	223	while(track) {
9f52c61f	224	fHistTracksPt[fCentBin]->Fill(track->Pt());
4bb8c6df	225	track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
9f52c61f	226	}
9f52c61f	227	}
46b40b53	228
8612dfc8	229	if (fCaloClustersCont) {
0925db21	230	AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
8612dfc8	231	while(cluster) {
9f52c61f	232	TLorentzVector nPart;
	233	cluster->GetMomentum(nPart, fVertex);
	234	fHistClustersPt[fCentBin]->Fill(nPart.Pt());
46b40b53	235	Double_t dx = cluster->GetTrackDx();
	236	Double_t dz = cluster->GetTrackDz();
	237	fHistClustDx->Fill(dx);
	238	fHistClustDz->Fill(dz);
0925db21	239	cluster = fCaloClustersCont->GetNextAcceptCluster();
9f52c61f	240	}
	241	}
	242
8612dfc8	243	if (fJetsCont) {
	244	AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
	245	while(jet) {
9239b066	246
9f52c61f	247	fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
	248	fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
	249
22339f55	250	Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
9f52c61f	251	fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
9f52c61f	252
388f9cd5	253	if (fHistJetsCorrPtArea[fCentBin]) {
	254	Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
	255	fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
	256	}
8612dfc8	257	jet = fJetsCont->GetNextAcceptJet();
9f52c61f	258	}
8612dfc8	259
8612dfc8	260	jet = fJetsCont->GetLeadingJet();
6b840003	261	if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
9f52c61f	262	}
9f52c61f	263
0925db21	264	CheckClusTrackMatching();
0925db21	265
9f52c61f	266	return kTRUE;
	267	}
	268
0925db21	269	//________________________________________________________________________
	270	void AliAnalysisTaskEmcalJetSample::CheckClusTrackMatching()
	271	{
	272
	273	if(!fTracksCont \|\| !fCaloClustersCont)
	274	return;
	275
	276	Double_t deta = 999;
	277	Double_t dphi = 999;
	278
	279	//Get closest cluster to track
	280	AliVTrack track = static_cast<AliVTrack>(fTracksCont->GetNextAcceptParticle(0));
	281	while(track) {
	282	//Get matched cluster
	283	Int_t emc1 = track->GetEMCALcluster();
	284	if(fCaloClustersCont && emc1>=0) {
	285	AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
	286	if(clusMatch) {
	287	AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
	288	fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
	289	}
	290	}
	291	track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
	292	}
	293
	294	//Get closest track to cluster
	295	AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
	296	while(cluster) {
	297	TLorentzVector nPart;
	298	cluster->GetMomentum(nPart, fVertex);
	299	fHistClustersPt[fCentBin]->Fill(nPart.Pt());
	300
	301	//Get matched track
	302	AliVTrack *mt = NULL;
	303	AliAODCaloCluster acl = dynamic_cast<AliAODCaloCluster>(cluster);
	304	if(acl) {
	305	if(acl->GetNTracksMatched()>1)
	306	mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
	307	}
	308	else {
	309	AliESDCaloCluster ecl = dynamic_cast<AliESDCaloCluster>(cluster);
	310	Int_t im = ecl->GetTrackMatchedIndex();
	311	if(fTracksCont && im>=0) {
	312	mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
	313	}
	314	}
	315	if(mt) {
	316	AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
	317	fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
	318
	319	/* //debugging
	320	if(mt->IsEMCAL()) {
	321	Int_t emc1 = mt->GetEMCALcluster();
	322	Printf("current id: %d emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
	323	AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
	324	AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
	325	Printf("deta: %f dphi: %f",deta,dphi);
	326	}
	327	*/
	328	}
	329	cluster = fCaloClustersCont->GetNextAcceptCluster();
	330	}
	331	}
	332
8612dfc8	333	//________________________________________________________________________
	334	void AliAnalysisTaskEmcalJetSample::ExecOnce() {
	335
	336	AliAnalysisTaskEmcalJet::ExecOnce();
	337
	338	if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
	339	if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
	340	if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
	341
	342	}
	343
9f52c61f	344	//________________________________________________________________________
	345	Bool_t AliAnalysisTaskEmcalJetSample::Run()
	346	{
	347	// Run analysis code here, if needed. It will be executed before FillHistograms().
	348
	349	return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
	350	}
	351
	352	//________________________________________________________________________
	353	void AliAnalysisTaskEmcalJetSample::Terminate(Option_t *)
	354	{
	355	// Called once at the end of the analysis.
	356	}