[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++) {
    fHistLeadingJetPhiEta[i] = 0;
    fHistLeadingJetPtArea[i] = 0;
    fHistLeadingJetCorrPtArea[i] = 0;
    fHistLeadingJetMCPtArea[i] = 0;
    fHistRhoVSleadJetPt[i] = 0;
    fHistJetPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
    fHistJetPtvsJetCorrPt[i] = 0;
    fHistJetsMCPtArea[i] = 0;
    fHistJetPtvsJetMCPt[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++) {
    fHistLeadingJetPhiEta[i] = 0;
    fHistLeadingJetPtArea[i] = 0;
    fHistLeadingJetCorrPtArea[i] = 0;
    fHistLeadingJetMCPtArea[i] = 0;
    fHistRhoVSleadJetPt[i] = 0;
    fHistJetPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
    fHistJetPtvsJetCorrPt[i] = 0;
    fHistJetsMCPtArea[i] = 0;
    fHistJetPtvsJetMCPt[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("No. 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(50, 0, 2);
  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 = "fHistLeadingJetPhiEta_";
    histname += i;
    fHistLeadingJetPhiEta[i] = new TH3F(histname.Data(), histname.Data(), 
					50, binsEta, 
					101, binsPhi, 
					nbinsZ, binsZ);
    fHistLeadingJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
    fHistLeadingJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
    fHistLeadingJetPhiEta[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
    fOutput->Add(fHistLeadingJetPhiEta[i]);

    histname = "fHistLeadingJetPtArea_";
    histname += i;
    fHistLeadingJetPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
					fNbins, binsPt, 
					50, binsArea,
					nbinsZ, binsZ);
    fHistLeadingJetPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
    fHistLeadingJetPtArea[i]->GetYaxis()->SetTitle("area");
    fHistLeadingJetPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
    fOutput->Add(fHistLeadingJetPtArea[i]);

    if (fIsEmbedded) {
      histname = "fHistLeadingJetMCPtArea_";
      histname += i;
      fHistLeadingJetMCPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
					    fNbins, binsPt, 
					    50, binsArea,
					    nbinsZ, binsZ);
      fHistLeadingJetMCPtArea[i]->GetXaxis()->SetTitle("p_{T,MC} (GeV/c)");
      fHistLeadingJetMCPtArea[i]->GetYaxis()->SetTitle("area");
      fHistLeadingJetMCPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
      fOutput->Add(fHistLeadingJetMCPtArea[i]);
    }

    if (!fRhoName.IsNull()) {
      histname = "fHistLeadingJetCorrPtArea_";
      histname += i;
      fHistLeadingJetCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
					      fNbins * 2, binsCorrPt, 
					      50, binsArea,
					      nbinsZ, binsZ);
      fHistLeadingJetCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
      fHistLeadingJetCorrPtArea[i]->GetYaxis()->SetTitle("area");
      fHistLeadingJetCorrPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
      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, 
				  50, 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]);

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

      histname = "fHistJetPtvsJetCorrPt_";
      histname += i;
      fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
					  fNbins, fMinBinPt, fMaxBinPt, 
					  fNbins * 2, -fMaxBinPt, fMaxBinPt);
      fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
      fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
      fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
      fOutput->Add(fHistJetPtvsJetCorrPt[i]);
    }

    if (fIsEmbedded) {
      histname = "fHistJetsMCPtArea_";
      histname += i;
      fHistJetsMCPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
				      fNbins, binsPt, 
				      50, binsArea,
				      nbinsZ, binsZ);
      fHistJetsMCPtArea[i]->GetXaxis()->SetTitle("p_{T,MC} (GeV/c)");
      fHistJetsMCPtArea[i]->GetYaxis()->SetTitle("area");
      fHistJetsMCPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
      fOutput->Add(fHistJetsMCPtArea[i]);

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

    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
    return kTRUE;

  static Int_t sortedJets[9999] = {-1};
  Bool_t r = GetSortedArray(sortedJets, fJets, fRhoVal);
  
  if (!r) // no accepted jets, skipping
    return kTRUE;

  // OK, event accepted

  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;

    Float_t ptLeading = GetLeadingHadronPt(jet);

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

    if (fIsEmbedded) 
      fHistLeadingJetMCPtArea[fCentBin]->Fill(jet->MCPt(), jet->Area(), ptLeading);      

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

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

    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);
    if (fHistJetPtvsJetCorrPt[fCentBin])
      fHistJetPtvsJetCorrPt[fCentBin]->Fill(jet->Pt(), corrPt);
    fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt(), ptLeading);

    if (fIsEmbedded) {
      fHistJetsMCPtArea[fCentBin]->Fill(jet->MCPt(), jet->Area(), ptLeading);
      fHistJetPtvsJetMCPt[fCentBin]->Fill(jet->Pt(), jet->MCPt());
    }

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