[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),
  fLeadingHadronType(0),
  fHistRhoVSleadJetPt(0),
  fNjetsVsCent(0)

{
  // Default constructor.

  for (Int_t i = 0; i < 4; i++) {
    fHistEvents[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHist2LeadingJetPt[i] = 0;
    fHistLeadingJetCorrPt[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),
  fLeadingHadronType(0),
  fHistRhoVSleadJetPt(0),
  fNjetsVsCent(0)
{
  // Standard constructor.

  for (Int_t i = 0; i < 4; i++) {
    fHistEvents[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHist2LeadingJetPt[i] = 0;
    fHistLeadingJetCorrPt[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.
}

//________________________________________________________________________
Float_t* AliAnalysisTaskSAJF::GenerateFixedBinArray(Int_t n, Float_t min, Float_t max) const
{
  Float_t *bins = new Float_t[n+1];

  Float_t binWidth = (max-min)/n;
  bins[0] = min;
  for (Int_t i = 1; i <= n; i++) {
    bins[i] = bins[i-1]+binWidth;
  }

  return bins;
}

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

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  fHistRhoVSleadJetPt = new TH2F("fHistRhoVSleadJetPt","fHistRhoVSleadJetPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
  fHistRhoVSleadJetPt->GetXaxis()->SetTitle("#rho * area (GeV/c)");
  fHistRhoVSleadJetPt->GetYaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
  fOutput->Add(fHistRhoVSleadJetPt);

  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 = "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->Fill(fRhoVal * jet->Area(), 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 = 0;

    if (fLeadingHadronType == 0)       // charged leading hadron
      ptLeading = jet->MaxTrackPt();
    else if (fLeadingHadronType == 1)  // neutral leading hadron
      ptLeading = jet->MaxClusterPt();
    else                               // charged or neutral
      ptLeading = jet->MaxPartPt();

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