[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskSAJF.cxx

// $Id$
//
// Jet 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 "AliVParticle.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"

#include "AliAnalysisTaskSAJF.h"

ClassImp(AliAnalysisTaskSAJF)

//________________________________________________________________________
AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF", kTRUE),
  fNjetsVsCent(0)

{
  // Default constructor.

  for (Int_t i = 0; i < 4; i++) {
    fHistEvents[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHist2LeadingJetPt[i] = 0;
    fHistLeadingJetCorrPt[i] = 0;
    fHistRhoVSleadJetPt[i] = 0;
    fHistJetPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
    fHistJetsNEFvsPt[i] = 0;
    fHistJetsZvsPt[i] = 0;
    fHistConstituents[i] = 0;
    fHistTracksJetPt[i] = 0;
    fHistClustersJetPt[i] = 0;
    fHistJetNconstVsPt[i] = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),
  fNjetsVsCent(0)
{
  // Standard constructor.

  for (Int_t i = 0; i < 4; i++) {
    fHistEvents[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHist2LeadingJetPt[i] = 0;
    fHistLeadingJetCorrPt[i] = 0;
    fHistRhoVSleadJetPt[i] = 0;
    fHistJetPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
    fHistJetsNEFvsPt[i] = 0;
    fHistJetsZvsPt[i] = 0;
    fHistConstituents[i] = 0;
    fHistTracksJetPt[i] = 0;
    fHistClustersJetPt[i] = 0;
    fHistJetNconstVsPt[i] = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

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

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

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  fNjetsVsCent = new TH2F("fNjetsVsCent","fNjetsVsCent", 100, 0, 100, 150, -0.5, 149.5);
  fNjetsVsCent->GetXaxis()->SetTitle("Centrality (%)");
  fNjetsVsCent->GetYaxis()->SetTitle("# of jets");
  fOutput->Add(fNjetsVsCent);

  const Int_t nbinsZ = 12;
  Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};

  Float_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
  Float_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
  Float_t *binsArea     = GenerateFixedBinArray(30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
  Float_t *binsEta      = GenerateFixedBinArray(50,-1, 1);
  Float_t *binsPhi      = GenerateFixedBinArray(101, 0, TMath::Pi() * 2.02);
  Float_t *bins120      = GenerateFixedBinArray(120, 0, 1.2);
  Float_t *bins150      = GenerateFixedBinArray(150, -0.5, 149.5);

  TString histname;

  for (Int_t i = 0; i < 4; i++) {
    histname = "fHistEvents_";
    histname += i;
    fHistEvents[i] = new TH1F(histname,histname, 6, 0, 6);
    fHistEvents[i]->GetXaxis()->SetTitle("Event state");
    fHistEvents[i]->GetYaxis()->SetTitle("counts");
    fHistEvents[i]->GetXaxis()->SetBinLabel(1, "No jets");
    fHistEvents[i]->GetXaxis()->SetBinLabel(2, "Max jet not found");
    fHistEvents[i]->GetXaxis()->SetBinLabel(3, "Rho == 0");
    fHistEvents[i]->GetXaxis()->SetBinLabel(4, "Max jet <= 0");
    fHistEvents[i]->GetXaxis()->SetBinLabel(5, "OK");
    fOutput->Add(fHistEvents[i]);

    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 = "fHist2LeadingJetPt_";
    histname += i;
    fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
    fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
    fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHist2LeadingJetPt[i]);

    histname = "fHistLeadingJetCorrPt_";
    histname += i;
    fHistLeadingJetCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
    fHistLeadingJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
    fHistLeadingJetCorrPt[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistLeadingJetCorrPt[i]);

    histname = "fHistRhoVSleadJetPt_";
    histname += i;
    fHistRhoVSleadJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt*2, fNbins, fMinBinPt, fMaxBinPt);
    fHistRhoVSleadJetPt[i]->GetXaxis()->SetTitle("#rho * area (GeV/c)");
    fHistRhoVSleadJetPt[i]->GetYaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
    fOutput->Add(fHistRhoVSleadJetPt[i]);
    
    histname = "fHistJetPhiEta_";
    histname += i;
    fHistJetPhiEta[i] = new TH3F(histname.Data(), histname.Data(), 
				 50, binsEta, 
				 101, binsPhi, 
				 nbinsZ, binsZ);
    fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
    fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
    fHistJetPhiEta[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
    fOutput->Add(fHistJetPhiEta[i]);
    
    histname = "fHistJetsPtArea_";
    histname += i;
    fHistJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
				  fNbins, binsPt, 
				  30, binsArea,
				  nbinsZ, binsZ);
    fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
    fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
    fHistJetsPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
    fOutput->Add(fHistJetsPtArea[i]);

    histname = "fHistJetsZvsPt_";
    histname += i;
    fHistJetsZvsPt[i] = new TH3F(histname.Data(), histname.Data(), 
				 120, bins120, 
				 fNbins, binsPt,
				 nbinsZ, binsZ);
    fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
    fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
    fHistJetsZvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
    fOutput->Add(fHistJetsZvsPt[i]);

    if (!fCaloName.IsNull()) {
      histname = "fHistJetsNEFvsPt_";
      histname += i;
      fHistJetsNEFvsPt[i] = new TH3F(histname.Data(), histname.Data(), 
				     120, bins120, 
				     fNbins, binsPt,
				     nbinsZ, binsZ);
      fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
      fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
      fHistJetsNEFvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
      fOutput->Add(fHistJetsNEFvsPt[i]);
    }

    histname = "fHistJetsCorrPtArea_";
    histname += i;
    fHistJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
				      fNbins * 2, binsCorrPt, 
				      30, binsArea,
				      nbinsZ, binsZ);
    fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
    fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
    fOutput->Add(fHistJetsCorrPtArea[i]);

    histname = "fHistConstituents_";
    histname += i;
    fHistConstituents[i] = new TH2F(histname.Data(), histname.Data(), 100, 1, 101, 100, -0.5, 99.5);
    fHistConstituents[i]->GetXaxis()->SetTitle("p_{T,part} (GeV/c)");
    fHistConstituents[i]->GetYaxis()->SetTitle("no. of particles");
    fHistConstituents[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistConstituents[i]);

    histname = "fHistTracksJetPt_";
    histname += i;
    fHistTracksJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
    fHistTracksJetPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
    fHistTracksJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
    fHistTracksJetPt[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistTracksJetPt[i]);

    if (!fCaloName.IsNull()) {
      histname = "fHistClustersJetPt_";
      histname += i;
      fHistClustersJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
      fHistClustersJetPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
      fHistClustersJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
      fHistClustersJetPt[i]->GetZaxis()->SetTitle("counts");
      fOutput->Add(fHistClustersJetPt[i]);
    }

    histname = "fHistJetNconstVsPt_";
    histname += i;
    fHistJetNconstVsPt[i] = new TH3F(histname.Data(), histname.Data(), 150, bins150, fNbins, binsPt, nbinsZ, binsZ);
    fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("# of constituents");
    fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
    fHistJetNconstVsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
    fOutput->Add(fHistJetNconstVsPt[i]);
  }

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

  delete binsPt;
  delete binsCorrPt;
  delete binsArea;
  delete binsEta;
  delete binsPhi;
  delete bins120;
}

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

  if (!fJets) {
    AliError(Form("%s - Jet array not provided, returning...", GetName()));
    return kFALSE;
  }

  if (fJets->GetEntriesFast() < 1) { // no jets in array, skipping
    fHistEvents[fCentBin]->Fill("No jets", 1);
    return kTRUE;
  }

  Int_t *sortedJets = GetSortedArray(fJets);
  
  if (!sortedJets || sortedJets[0] < 0) { // no accepted jets, skipping
    fHistEvents[fCentBin]->Fill("No jets", 1);
    return kTRUE;
  }

  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
  if (!jet) {  // error, I cannot get the leading jet from collection (should never happen), skipping
    fHistEvents[fCentBin]->Fill("Max jet not found", 1);
    return kTRUE;
  }

  // OK, event accepted

  Float_t maxJetCorrPt = jet->Pt() - fRhoVal * jet->Area();

  if (fRhoVal == 0) 
    fHistEvents[fCentBin]->Fill("Rho == 0", 1);

  else if (maxJetCorrPt <= 0)
    fHistEvents[fCentBin]->Fill("Max jet <= 0", 1);

  else
    fHistEvents[fCentBin]->Fill("OK", 1);

  if (jet) {
    fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
    fHistRhoVSleadJetPt[fCentBin]->Fill(fRhoVal, jet->Pt());
    fHistLeadingJetCorrPt[fCentBin]->Fill(maxJetCorrPt);
  }

  AliEmcalJet* jet2 = 0;
  if (sortedJets[1] >= 0)
    jet2 = static_cast<AliEmcalJet*>(fJets->At(sortedJets[1]));

  if (jet2)
    fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());

  Int_t njets = DoJetLoop();

  fNjetsVsCent->Fill(fCent, njets);

  return kTRUE;
}

//________________________________________________________________________
Int_t AliAnalysisTaskSAJF::DoJetLoop()
{
  // Do the jet loop.

  if (!fJets)
    return 0;

  const Int_t njets = fJets->GetEntriesFast();
  Int_t nAccJets = 0;

  TH1F constituents("constituents", "constituents", 
		    fHistConstituents[0]->GetNbinsX(), fHistConstituents[0]->GetXaxis()->GetXmin(), fHistConstituents[0]->GetXaxis()->GetXmax()); 

  for (Int_t ij = 0; ij < njets; ij++) {

    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));

    if (!jet) {
      AliError(Form("Could not receive jet %d", ij));
      continue;
    }  

    if (!AcceptJet(jet))
      continue;

    Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();

    Float_t ptLeading = GetLeadingHadronPt(jet);

    fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi(), ptLeading);
    fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area(), ptLeading);
    fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area(), ptLeading);
    fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt(), ptLeading);

    if (fCaloClusters)
      fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt(), ptLeading);

    if (fTracks) {
      for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
	AliVParticle *track = jet->TrackAt(it, fTracks);
	if (track) {
	  fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt(), ptLeading);
	  constituents.Fill(track->Pt());
	  fHistTracksJetPt[fCentBin]->Fill(track->Pt(), jet->Pt());
	}
      }
    }

    if (fCaloClusters) {
      for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
	AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
	
	if (cluster) {
	  TLorentzVector nPart;
	  cluster->GetMomentum(nPart, fVertex);
	  fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt(), ptLeading);
	  constituents.Fill(nPart.Pt());
	  fHistClustersJetPt[fCentBin]->Fill(nPart.Pt(), jet->Pt());
	}
      }
    }

    for (Int_t i = 1; i <= constituents.GetNbinsX(); i++) {
      fHistConstituents[fCentBin]->Fill(constituents.GetBinCenter(i), constituents.GetBinContent(i));
    }

    constituents.Reset();
    nAccJets++;
  } //jet loop 

  return nAccJets;
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
	1	// $Id$
	2	//
	3	// Jet analysis task.
	4	//
	5	// Author: S.Aiola
	6
	7	#include <TClonesArray.h>
	8	#include <TH1F.h>
	9	#include <TH2F.h>
	10	#include <TH3F.h>
	11	#include <TList.h>
	12	#include <TLorentzVector.h>
	13
	14	#include "AliVCluster.h"
	15	#include "AliVParticle.h"
	16	#include "AliVTrack.h"
	17	#include "AliEmcalJet.h"
	18	#include "AliRhoParameter.h"
	19	#include "AliLog.h"
	20
	21	#include "AliAnalysisTaskSAJF.h"
	22
	23	ClassImp(AliAnalysisTaskSAJF)
	24
	25	//________________________________________________________________________
	26	AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
	27	AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF", kTRUE),
	28	fNjetsVsCent(0)
	29
	30	{
	31	// Default constructor.
	32
	33	for (Int_t i = 0; i < 4; i++) {
	34	fHistEvents[i] = 0;
	35	fHistLeadingJetPt[i] = 0;
	36	fHist2LeadingJetPt[i] = 0;
	37	fHistLeadingJetCorrPt[i] = 0;
	38	fHistRhoVSleadJetPt[i] = 0;
	39	fHistJetPhiEta[i] = 0;
	40	fHistJetsPtArea[i] = 0;
	41	fHistJetsCorrPtArea[i] = 0;
	42	fHistJetsNEFvsPt[i] = 0;
	43	fHistJetsZvsPt[i] = 0;
	44	fHistConstituents[i] = 0;
	45	fHistTracksJetPt[i] = 0;
	46	fHistClustersJetPt[i] = 0;
	47	fHistJetNconstVsPt[i] = 0;
	48	}
	49
	50	SetMakeGeneralHistograms(kTRUE);
	51	}
	52
	53	//________________________________________________________________________
	54	AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
	55	AliAnalysisTaskEmcalJet(name, kTRUE),
	56	fNjetsVsCent(0)
	57	{
	58	// Standard constructor.
	59
	60	for (Int_t i = 0; i < 4; i++) {
	61	fHistEvents[i] = 0;
	62	fHistLeadingJetPt[i] = 0;
	63	fHist2LeadingJetPt[i] = 0;
	64	fHistLeadingJetCorrPt[i] = 0;
	65	fHistRhoVSleadJetPt[i] = 0;
	66	fHistJetPhiEta[i] = 0;
	67	fHistJetsPtArea[i] = 0;
	68	fHistJetsCorrPtArea[i] = 0;
	69	fHistJetsNEFvsPt[i] = 0;
	70	fHistJetsZvsPt[i] = 0;
	71	fHistConstituents[i] = 0;
	72	fHistTracksJetPt[i] = 0;
	73	fHistClustersJetPt[i] = 0;
	74	fHistJetNconstVsPt[i] = 0;
	75	}
	76
	77	SetMakeGeneralHistograms(kTRUE);
	78	}
	79
	80	//________________________________________________________________________
	81	AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
	82	{
	83	// Destructor.
	84	}
	85
	86	//________________________________________________________________________
	87	void AliAnalysisTaskSAJF::UserCreateOutputObjects()
	88	{
	89	// Create user output.
	90
	91	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
	92
	93	fNjetsVsCent = new TH2F("fNjetsVsCent","fNjetsVsCent", 100, 0, 100, 150, -0.5, 149.5);
	94	fNjetsVsCent->GetXaxis()->SetTitle("Centrality (%)");
	95	fNjetsVsCent->GetYaxis()->SetTitle("# of jets");
	96	fOutput->Add(fNjetsVsCent);
	97
	98	const Int_t nbinsZ = 12;
	99	Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
	100
	101	Float_t *binsPt = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
	102	Float_t binsCorrPt = GenerateFixedBinArray(fNbins2, -fMaxBinPt, fMaxBinPt);
	103	Float_t binsArea = GenerateFixedBinArray(30, 0, fJetRadius fJetRadius * TMath::Pi() * 3);
	104	Float_t *binsEta = GenerateFixedBinArray(50,-1, 1);
	105	Float_t binsPhi = GenerateFixedBinArray(101, 0, TMath::Pi() 2.02);
	106	Float_t *bins120 = GenerateFixedBinArray(120, 0, 1.2);
	107	Float_t *bins150 = GenerateFixedBinArray(150, -0.5, 149.5);
	108
	109	TString histname;
	110
	111	for (Int_t i = 0; i < 4; i++) {
	112	histname = "fHistEvents_";
	113	histname += i;
	114	fHistEvents[i] = new TH1F(histname,histname, 6, 0, 6);
	115	fHistEvents[i]->GetXaxis()->SetTitle("Event state");
	116	fHistEvents[i]->GetYaxis()->SetTitle("counts");
	117	fHistEvents[i]->GetXaxis()->SetBinLabel(1, "No jets");
	118	fHistEvents[i]->GetXaxis()->SetBinLabel(2, "Max jet not found");
	119	fHistEvents[i]->GetXaxis()->SetBinLabel(3, "Rho == 0");
	120	fHistEvents[i]->GetXaxis()->SetBinLabel(4, "Max jet <= 0");
	121	fHistEvents[i]->GetXaxis()->SetBinLabel(5, "OK");
	122	fOutput->Add(fHistEvents[i]);
	123
	124	histname = "fHistLeadingJetPt_";
	125	histname += i;
	126	fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
	127	fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
	128	fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
	129	fOutput->Add(fHistLeadingJetPt[i]);
	130
	131	histname = "fHist2LeadingJetPt_";
	132	histname += i;
	133	fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
	134	fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
	135	fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
	136	fOutput->Add(fHist2LeadingJetPt[i]);
	137
	138	histname = "fHistLeadingJetCorrPt_";
	139	histname += i;
	140	fHistLeadingJetCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
	141	fHistLeadingJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
	142	fHistLeadingJetCorrPt[i]->GetYaxis()->SetTitle("counts");
	143	fOutput->Add(fHistLeadingJetCorrPt[i]);
	144
	145	histname = "fHistRhoVSleadJetPt_";
	146	histname += i;
	147	fHistRhoVSleadJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt*2, fNbins, fMinBinPt, fMaxBinPt);
	148	fHistRhoVSleadJetPt[i]->GetXaxis()->SetTitle("#rho * area (GeV/c)");
	149	fHistRhoVSleadJetPt[i]->GetYaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
	150	fOutput->Add(fHistRhoVSleadJetPt[i]);
	151
	152	histname = "fHistJetPhiEta_";
	153	histname += i;
	154	fHistJetPhiEta[i] = new TH3F(histname.Data(), histname.Data(),
	155	50, binsEta,
	156	101, binsPhi,
	157	nbinsZ, binsZ);
	158	fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
	159	fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
	160	fHistJetPhiEta[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
	161	fOutput->Add(fHistJetPhiEta[i]);
	162
	163	histname = "fHistJetsPtArea_";
	164	histname += i;
	165	fHistJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(),
	166	fNbins, binsPt,
	167	30, binsArea,
	168	nbinsZ, binsZ);
	169	fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
	170	fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
	171	fHistJetsPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
	172	fOutput->Add(fHistJetsPtArea[i]);
	173
	174	histname = "fHistJetsZvsPt_";
	175	histname += i;
	176	fHistJetsZvsPt[i] = new TH3F(histname.Data(), histname.Data(),
	177	120, bins120,
	178	fNbins, binsPt,
	179	nbinsZ, binsZ);
	180	fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
	181	fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
	182	fHistJetsZvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
	183	fOutput->Add(fHistJetsZvsPt[i]);
	184
	185	if (!fCaloName.IsNull()) {
	186	histname = "fHistJetsNEFvsPt_";
	187	histname += i;
	188	fHistJetsNEFvsPt[i] = new TH3F(histname.Data(), histname.Data(),
	189	120, bins120,
	190	fNbins, binsPt,
	191	nbinsZ, binsZ);
	192	fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
	193	fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
	194	fHistJetsNEFvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
	195	fOutput->Add(fHistJetsNEFvsPt[i]);
	196	}
	197
	198	histname = "fHistJetsCorrPtArea_";
	199	histname += i;
	200	fHistJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(),
	201	fNbins * 2, binsCorrPt,
	202	30, binsArea,
	203	nbinsZ, binsZ);
	204	fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
	205	fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
	206	fOutput->Add(fHistJetsCorrPtArea[i]);
	207
	208	histname = "fHistConstituents_";
	209	histname += i;
	210	fHistConstituents[i] = new TH2F(histname.Data(), histname.Data(), 100, 1, 101, 100, -0.5, 99.5);
	211	fHistConstituents[i]->GetXaxis()->SetTitle("p_{T,part} (GeV/c)");
	212	fHistConstituents[i]->GetYaxis()->SetTitle("no. of particles");
	213	fHistConstituents[i]->GetZaxis()->SetTitle("counts");
	214	fOutput->Add(fHistConstituents[i]);
	215
	216	histname = "fHistTracksJetPt_";
	217	histname += i;
	218	fHistTracksJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
	219	fHistTracksJetPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
	220	fHistTracksJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
	221	fHistTracksJetPt[i]->GetZaxis()->SetTitle("counts");
	222	fOutput->Add(fHistTracksJetPt[i]);
	223
	224	if (!fCaloName.IsNull()) {
	225	histname = "fHistClustersJetPt_";
	226	histname += i;
	227	fHistClustersJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
	228	fHistClustersJetPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
	229	fHistClustersJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
	230	fHistClustersJetPt[i]->GetZaxis()->SetTitle("counts");
	231	fOutput->Add(fHistClustersJetPt[i]);
	232	}
	233
	234	histname = "fHistJetNconstVsPt_";
	235	histname += i;
	236	fHistJetNconstVsPt[i] = new TH3F(histname.Data(), histname.Data(), 150, bins150, fNbins, binsPt, nbinsZ, binsZ);
	237	fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("# of constituents");
	238	fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
	239	fHistJetNconstVsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
	240	fOutput->Add(fHistJetNconstVsPt[i]);
	241	}
	242
	243	PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
	244
	245	delete binsPt;
	246	delete binsCorrPt;
	247	delete binsArea;
	248	delete binsEta;
	249	delete binsPhi;
	250	delete bins120;
	251	}
	252
	253	//________________________________________________________________________
	254	Bool_t AliAnalysisTaskSAJF::FillHistograms()
	255	{
	256	// Fill histograms.
	257
	258	if (!fJets) {
	259	AliError(Form("%s - Jet array not provided, returning...", GetName()));
	260	return kFALSE;
	261	}
	262
	263	if (fJets->GetEntriesFast() < 1) { // no jets in array, skipping
	264	fHistEvents[fCentBin]->Fill("No jets", 1);
	265	return kTRUE;
	266	}
	267
	268	Int_t *sortedJets = GetSortedArray(fJets);
	269
	270	if (!sortedJets \|\| sortedJets[0] < 0) { // no accepted jets, skipping
	271	fHistEvents[fCentBin]->Fill("No jets", 1);
	272	return kTRUE;
	273	}
	274
	275	AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
	276	if (!jet) { // error, I cannot get the leading jet from collection (should never happen), skipping
	277	fHistEvents[fCentBin]->Fill("Max jet not found", 1);
	278	return kTRUE;
	279	}
	280
	281	// OK, event accepted
	282
	283	Float_t maxJetCorrPt = jet->Pt() - fRhoVal * jet->Area();
	284
	285	if (fRhoVal == 0)
	286	fHistEvents[fCentBin]->Fill("Rho == 0", 1);
	287
	288	else if (maxJetCorrPt <= 0)
	289	fHistEvents[fCentBin]->Fill("Max jet <= 0", 1);
	290
	291	else
	292	fHistEvents[fCentBin]->Fill("OK", 1);
	293
	294	if (jet) {
	295	fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
	296	fHistRhoVSleadJetPt[fCentBin]->Fill(fRhoVal, jet->Pt());
	297	fHistLeadingJetCorrPt[fCentBin]->Fill(maxJetCorrPt);
	298	}
	299
	300	AliEmcalJet* jet2 = 0;
	301	if (sortedJets[1] >= 0)
	302	jet2 = static_cast<AliEmcalJet*>(fJets->At(sortedJets[1]));
	303
	304	if (jet2)
	305	fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
	306
	307	Int_t njets = DoJetLoop();
	308
	309	fNjetsVsCent->Fill(fCent, njets);
	310
	311	return kTRUE;
	312	}
	313
	314	//________________________________________________________________________
	315	Int_t AliAnalysisTaskSAJF::DoJetLoop()
	316	{
	317	// Do the jet loop.
	318
	319	if (!fJets)
	320	return 0;
	321
	322	const Int_t njets = fJets->GetEntriesFast();
	323	Int_t nAccJets = 0;
	324
	325	TH1F constituents("constituents", "constituents",
	326	fHistConstituents[0]->GetNbinsX(), fHistConstituents[0]->GetXaxis()->GetXmin(), fHistConstituents[0]->GetXaxis()->GetXmax());
	327
	328	for (Int_t ij = 0; ij < njets; ij++) {
	329
	330	AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
	331
	332	if (!jet) {
	333	AliError(Form("Could not receive jet %d", ij));
	334	continue;
	335	}
	336
	337	if (!AcceptJet(jet))
	338	continue;
	339
	340	Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
	341
	342	Float_t ptLeading = GetLeadingHadronPt(jet);
	343
	344	fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi(), ptLeading);
	345	fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area(), ptLeading);
	346	fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area(), ptLeading);
	347	fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt(), ptLeading);
	348
	349	if (fCaloClusters)
	350	fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt(), ptLeading);
	351
	352	if (fTracks) {
	353	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
	354	AliVParticle *track = jet->TrackAt(it, fTracks);
	355	if (track) {
	356	fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt(), ptLeading);
	357	constituents.Fill(track->Pt());
	358	fHistTracksJetPt[fCentBin]->Fill(track->Pt(), jet->Pt());
	359	}
	360	}
	361	}
	362
	363	if (fCaloClusters) {
	364	for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
	365	AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
	366
	367	if (cluster) {
	368	TLorentzVector nPart;
	369	cluster->GetMomentum(nPart, fVertex);
	370	fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt(), ptLeading);
	371	constituents.Fill(nPart.Pt());
	372	fHistClustersJetPt[fCentBin]->Fill(nPart.Pt(), jet->Pt());
	373	}
	374	}
	375	}
	376
	377	for (Int_t i = 1; i <= constituents.GetNbinsX(); i++) {
	378	fHistConstituents[fCentBin]->Fill(constituents.GetBinCenter(i), constituents.GetBinContent(i));
	379	}
	380
	381	constituents.Reset();
	382	nAccJets++;
	383	} //jet loop
	384
	385	return nAccJets;
	386	}
	387
	388	//________________________________________________________________________
	389	void AliAnalysisTaskSAJF::Terminate(Option_t *)
	390	{
	391	// Called once at the end of the analysis.
	392	}