[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;
    fHistLeadingJetPhiEta[i] = 0;
    fHistLeadingJetPtArea[i] = 0;
    fHistLeadingJetCorrPtArea[i] = 0;
    fHistRhoVSleadJetPt[i] = 0;
    fHistJetPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
    fHistJetsNEFvsPt[i] = 0;
    fHistJetsCEFvsCEFPt[i] = 0;
    fHistJetsZvsPt[i] = 0;
    fHistConstituents[i] = 0;
    fHistTracksJetPt[i] = 0;
    fHistClustersJetPt[i] = 0;
    fHistTracksPtDist[i] = 0;
    fHistClustersPtDist[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;
    fHistLeadingJetPhiEta[i] = 0;
    fHistLeadingJetPtArea[i] = 0;
    fHistLeadingJetCorrPtArea[i] = 0;
    fHistRhoVSleadJetPt[i] = 0;
    fHistJetPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
    fHistJetsNEFvsPt[i] = 0;
    fHistJetsCEFvsCEFPt[i] = 0;
    fHistJetsZvsPt[i] = 0;
    fHistConstituents[i] = 0;
    fHistTracksJetPt[i] = 0;
    fHistClustersJetPt[i] = 0;
    fHistTracksPtDist[i] = 0;
    fHistClustersPtDist[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 = "fHistLeadingJetPhiEta_";
    histname += i;
    fHistLeadingJetPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50,-1, 1, 101, 0, TMath::Pi() * 2.02);
    fHistLeadingJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
    fHistLeadingJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
    fOutput->Add(fHistLeadingJetPhiEta[i]);

    histname = "fHistLeadingJetPtArea_";
    histname += i;
    fHistLeadingJetPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
    fHistLeadingJetPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
    fHistLeadingJetPtArea[i]->GetYaxis()->SetTitle("area");
    fOutput->Add(fHistLeadingJetPtArea[i]);

    if (!fRhoName.IsNull()) {
      histname = "fHistLeadingJetCorrPtArea_";
      histname += i;
      fHistLeadingJetCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, 30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
      fHistLeadingJetCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
      fHistLeadingJetCorrPtArea[i]->GetYaxis()->SetTitle("area");
      fOutput->Add(fHistLeadingJetCorrPtArea[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]);

    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 = "fHistJetsCEFvsCEFPt_";
    histname += i;
    fHistJetsCEFvsCEFPt[i] = new TH3F(histname.Data(), histname.Data(), 
				      120, bins120, 
				      fNbins, binsPt,
				      nbinsZ, binsZ);
    fHistJetsCEFvsCEFPt[i]->GetXaxis()->SetTitle("1-NEF");
    fHistJetsCEFvsCEFPt[i]->GetYaxis()->SetTitle("(1-NEF)*p_{T}^{raw} (GeV/c)");
    fHistJetsCEFvsCEFPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
    fOutput->Add(fHistJetsCEFvsCEFPt[i]);

    if (!fRhoName.IsNull()) {
      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]);

    histname = "fHistTracksPtDist_";
    histname += i;
    fHistTracksPtDist[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, 100, 0, 5);
    fHistTracksPtDist[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
    fHistTracksPtDist[i]->GetYaxis()->SetTitle("d");
    fHistTracksPtDist[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistTracksPtDist[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 = "fHistClustersPtDist_";
      histname += i;
      fHistClustersPtDist[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, 100, 0, 5);
      fHistClustersPtDist[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
      fHistClustersPtDist[i]->GetYaxis()->SetTitle("d");
      fHistClustersPtDist[i]->GetZaxis()->SetTitle("counts");
      fOutput->Add(fHistClustersPtDist[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;
  }

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

  // OK, event accepted

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

  for (Int_t i = 0; i < fNLeadingJets && i < fJets->GetEntriesFast(); i++) {
    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[i]));

    if (!jet) {
      AliError(Form("Could not receive jet %d", sortedJets[i]));
      continue;
    }  

    if (!AcceptJet(jet))
      continue;

    fHistLeadingJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
    fHistLeadingJetPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());

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

    if (fHistLeadingJetCorrPtArea[fCentBin])
      fHistLeadingJetCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());

    if (i==0 && fHistRhoVSleadJetPt[fCentBin]) 
      fHistRhoVSleadJetPt[fCentBin]->Fill(fRhoVal, jet->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);
    if (fHistJetsCorrPtArea[fCentBin])
      fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area(), ptLeading);
    fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt(), ptLeading);

    fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt(), ptLeading);
    fHistJetsCEFvsCEFPt[fCentBin]->Fill(1-jet->NEF(), (1-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());
	  Double_t dist = TMath::Sqrt((track->Eta() - jet->Eta()) * (track->Eta() - jet->Eta()) + (track->Phi() - jet->Phi()) * (track->Phi() - jet->Phi()));
	  fHistTracksPtDist[fCentBin]->Fill(track->Pt(), dist);
	}
      }
    }

    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());
	  Double_t dist = TMath::Sqrt((nPart.Eta() - jet->Eta()) * (nPart.Eta() - jet->Eta()) + (nPart.Phi() - jet->Phi()) * (nPart.Phi() - jet->Phi()));
	  fHistClustersPtDist[fCentBin]->Fill(nPart.Pt(), dist);
	}
      }
    }

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