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

// $Id$
//
// Jet embedding deltaPt task.
//
// Author: M. Verweij

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

#include "AliVCluster.h"
#include "AliVParticle.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliJetContainer.h"
#include "AliParticleContainer.h"
#include "AliClusterContainer.h"

#include "AliAnalysisTaskDeltaPtJEmb.h"

ClassImp(AliAnalysisTaskDeltaPtJEmb)

//________________________________________________________________________
AliAnalysisTaskDeltaPtJEmb::AliAnalysisTaskDeltaPtJEmb() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPtJEmb", kTRUE),
  fJetsCont(0),
  fTracksCont(0),
  fCaloClustersCont(0),
  fMinFractionShared(0.5),
  fHistEmbJetsPtArea(0),
  fHistEmbJetsCorrPtArea(0),
  fHistEmbPartPtvsJetPt(0),
  fHistEmbPartPtvsJetCorrPt(0),
  fHistJetPtvsJetCorrPt(0),
  fHistDistLeadPart2JetAxis(0),
  fHistEmbBkgArea(0),
  fHistRhoVSEmbBkg(0),
  fHistDeltaPtEmbArea(0),
  fHistDeltaPtEmbvsEP(0),
  fHistDeltaPtEmbPtProbe(0),
  fHistEmbJetsPhiEta(0),
  fHistLeadPartPhiEta(0)
{
  // Default constructor.

  fHistEmbJetsPtArea = new TH3*[fNcentBins];
  fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
  fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
  fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
  fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
  fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
  fHistEmbBkgArea = new TH2*[fNcentBins];
  fHistRhoVSEmbBkg = new TH2*[fNcentBins];
  fHistDeltaPtEmbArea = new TH2*[fNcentBins];
  fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
  fHistDeltaPtEmbPtProbe = new TH2*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fHistEmbJetsPtArea[i] = 0;
    fHistEmbJetsCorrPtArea[i] = 0;
    fHistEmbPartPtvsJetPt[i] = 0;
    fHistEmbPartPtvsJetCorrPt[i] = 0;
    fHistJetPtvsJetCorrPt[i] = 0;
    fHistDistLeadPart2JetAxis[i] = 0;
    fHistEmbBkgArea[i] = 0;
    fHistRhoVSEmbBkg[i] = 0;
    fHistDeltaPtEmbArea[i] = 0;
    fHistDeltaPtEmbvsEP[i] = 0;
    fHistDeltaPtEmbPtProbe[i] = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskDeltaPtJEmb::AliAnalysisTaskDeltaPtJEmb(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),
  fJetsCont(0),
  fTracksCont(0),
  fCaloClustersCont(0),
  fMinFractionShared(0.5),
  fHistEmbJetsPtArea(0),
  fHistEmbJetsCorrPtArea(0),
  fHistEmbPartPtvsJetPt(0),
  fHistEmbPartPtvsJetCorrPt(0),
  fHistJetPtvsJetCorrPt(0),
  fHistDistLeadPart2JetAxis(0),
  fHistEmbBkgArea(0),
  fHistRhoVSEmbBkg(0),
  fHistDeltaPtEmbArea(0),
  fHistDeltaPtEmbvsEP(0),
  fHistDeltaPtEmbPtProbe(0),
  fHistEmbJetsPhiEta(0),
  fHistLeadPartPhiEta(0)
{
  // Standard constructor.

  fHistEmbJetsPtArea = new TH3*[fNcentBins];
  fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
  fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
  fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
  fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
  fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
  fHistEmbBkgArea = new TH2*[fNcentBins];
  fHistRhoVSEmbBkg = new TH2*[fNcentBins];
  fHistDeltaPtEmbArea = new TH2*[fNcentBins];
  fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
  fHistDeltaPtEmbPtProbe = new TH2*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fHistEmbJetsPtArea[i] = 0;
    fHistEmbJetsCorrPtArea[i] = 0;
    fHistEmbPartPtvsJetPt[i] = 0;
    fHistEmbPartPtvsJetCorrPt[i] = 0;
    fHistJetPtvsJetCorrPt[i] = 0;
    fHistDistLeadPart2JetAxis[i] = 0;
    fHistEmbBkgArea[i] = 0;
    fHistRhoVSEmbBkg[i] = 0;
    fHistDeltaPtEmbArea[i] = 0;
    fHistDeltaPtEmbvsEP[i] = 0;
    fHistDeltaPtEmbPtProbe[i] = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

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

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  fJetsCont = GetJetContainer("Jets");
  fTracksCont = GetParticleContainer("Tracks");
  fCaloClustersCont = GetClusterContainer("CaloClusters");
  
  fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
  fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
  fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistEmbJetsPhiEta);
  
  fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
  fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
  fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistLeadPartPhiEta);
  
  TString histname;

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

  Double_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
  Double_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
  Double_t *binsArea     = GenerateFixedBinArray(50, 0, 2);

  for (Int_t i = 0; i < fNcentBins; i++) {
    histname = "fHistEmbJetsPtArea_";
    histname += i;
    fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
    fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
    fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
    fOutput->Add(fHistEmbJetsPtArea[i]);

    histname = "fHistEmbJetsCorrPtArea_";
    histname += i;
    fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
    fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
    fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
    fOutput->Add(fHistEmbJetsCorrPtArea[i]);

    histname = "fHistEmbPartPtvsJetPt_";
    histname += i;
    fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
    fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
    fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
    fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistEmbPartPtvsJetPt[i]);

    histname = "fHistEmbPartPtvsJetCorrPt_";
    histname += i;
    fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), 
					    fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
    fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
    fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
    fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);

    histname = "fHistJetPtvsJetCorrPt_";
    histname += i;
    fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), 
					fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
    fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
    fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
    fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistJetPtvsJetCorrPt[i]);

    histname = "fHistDistLeadPart2JetAxis_";
    histname += i;
    fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
    fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
    fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistDistLeadPart2JetAxis[i]);

    histname = "fHistEmbBkgArea_";
    histname += i;
    fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
    fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
    fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
    fOutput->Add(fHistEmbBkgArea[i]);

    histname = "fHistRhoVSEmbBkg_";
    histname += i;
    fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
    fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
    fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
    fOutput->Add(fHistRhoVSEmbBkg[i]);
      
    histname = "fHistDeltaPtEmbArea_";
    histname += i;
    fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(), 
				      50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
    fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
    fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
    fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistDeltaPtEmbArea[i]);

    histname = "fHistDeltaPtEmbvsEP_";
    histname += i;
    fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(), 
				      402, -TMath::Pi()*1.01, TMath::Pi()*3.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
    fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}");
    fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
    fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistDeltaPtEmbvsEP[i]);

    histname = "fHistDeltaPtEmbPtProbe_";
    histname += i;
    fHistDeltaPtEmbPtProbe[i] = new TH2F(histname.Data(), histname.Data(), 
					 fNbins, fMinBinPt, fMaxBinPt, fNbins * 2, -fMaxBinPt, fMaxBinPt);
    fHistDeltaPtEmbPtProbe[i]->GetXaxis()->SetTitle("#it{p}_{T,probe} (GeV/#it{c})");
    fHistDeltaPtEmbPtProbe[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
    fHistDeltaPtEmbPtProbe[i]->GetZaxis()->SetTitle("counts");
    fOutput->Add(fHistDeltaPtEmbPtProbe[i]);
  }

  delete[] binsPt;
  delete[] binsCorrPt;
  delete[] binsArea;

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

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

  // ************
  // Embedding
  // Jets should already be matched, GetClosestJet gives matched partner. Cut on shared fraction to be applied
  // _________________________________

  if (fJetsCont) {
    AliEmcalJet *jet1 = NULL;
    AliEmcalJet *jet2 = NULL;
    fJetsCont->ResetCurrentID();
    while((jet1 = fJetsCont->GetNextAcceptJet())) {

      jet2 = jet1->ClosestJet();
      if(!jet2) continue;

      Double_t fraction = fJetsCont->GetFractionSharedPt(jet1);
      if(fMinFractionShared>0. && fraction<fMinFractionShared) continue;

      TLorentzVector mom;
      fJetsCont->GetLeadingHadronMomentum(mom,jet1);
      
      Double_t distLeading2Jet = TMath::Sqrt((jet1->Eta() - mom.Eta()) * (jet1->Eta() - mom.Eta()) + (jet1->Phi() - mom.Phi()+TMath::Pi()) * (jet1->Phi() - mom.Phi()+TMath::Pi()));
      
      fHistEmbPartPtvsJetPt[fCentBin]->Fill(jet2->Pt(), jet1->Pt());
      fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(jet2->Pt(), jet1->Pt() - jet1->Area() * fRhoVal);
      fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi()+TMath::Pi());
      fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
      
      fHistEmbJetsPtArea[fCentBin]->Fill(jet1->Area(), jet1->Pt(), mom.Pt());
      fHistEmbJetsCorrPtArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area(), mom.Pt());
      fHistEmbJetsPhiEta->Fill(jet1->Eta(), jet1->Phi());
      fHistJetPtvsJetCorrPt[fCentBin]->Fill(jet1->Pt(), jet1->Pt() - fRhoVal * jet1->Area());
      
      fHistEmbBkgArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - jet2->Pt());
      fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * jet1->Area(), jet1->Pt() - jet2->Pt());
      fHistDeltaPtEmbArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
      fHistDeltaPtEmbvsEP[fCentBin]->Fill(jet1->Phi() - fEPV0, jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
      fHistDeltaPtEmbPtProbe[fCentBin]->Fill(jet2->Pt(), jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
    }
  }

  return kTRUE;
}

//________________________________________________________________________
void AliAnalysisTaskDeltaPtJEmb::ExecOnce()
{
  // Initialize the analysis.

  AliAnalysisTaskEmcalJet::ExecOnce();

  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
}
Commit	Line	Data
3a02e8fb	1	// $Id$
	2	//
	3	// Jet embedding deltaPt task.
	4	//
	5	// Author: M. Verweij
	6
	7	#include <TClonesArray.h>
	8	#include <TH1F.h>
	9	#include <TH2F.h>
	10	#include <TH3F.h>
	11	#include <TList.h>
	12	#include <TLorentzVector.h>
	13	#include <TRandom3.h>
	14
	15	#include "AliVCluster.h"
	16	#include "AliVParticle.h"
	17	#include "AliVTrack.h"
	18	#include "AliEmcalJet.h"
	19	#include "AliRhoParameter.h"
	20	#include "AliLog.h"
	21	#include "AliJetContainer.h"
	22	#include "AliParticleContainer.h"
	23	#include "AliClusterContainer.h"
	24
	25	#include "AliAnalysisTaskDeltaPtJEmb.h"
	26
	27	ClassImp(AliAnalysisTaskDeltaPtJEmb)
	28
	29	//________________________________________________________________________
	30	AliAnalysisTaskDeltaPtJEmb::AliAnalysisTaskDeltaPtJEmb() :
	31	AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPtJEmb", kTRUE),
	32	fJetsCont(0),
	33	fTracksCont(0),
	34	fCaloClustersCont(0),
	35	fMinFractionShared(0.5),
	36	fHistEmbJetsPtArea(0),
	37	fHistEmbJetsCorrPtArea(0),
	38	fHistEmbPartPtvsJetPt(0),
	39	fHistEmbPartPtvsJetCorrPt(0),
	40	fHistJetPtvsJetCorrPt(0),
	41	fHistDistLeadPart2JetAxis(0),
	42	fHistEmbBkgArea(0),
	43	fHistRhoVSEmbBkg(0),
	44	fHistDeltaPtEmbArea(0),
	45	fHistDeltaPtEmbvsEP(0),
	46	fHistDeltaPtEmbPtProbe(0),
	47	fHistEmbJetsPhiEta(0),
	48	fHistLeadPartPhiEta(0)
	49	{
	50	// Default constructor.
	51
	52	fHistEmbJetsPtArea = new TH3*[fNcentBins];
	53	fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
	54	fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
	55	fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
	56	fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
	57	fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
	58	fHistEmbBkgArea = new TH2*[fNcentBins];
	59	fHistRhoVSEmbBkg = new TH2*[fNcentBins];
	60	fHistDeltaPtEmbArea = new TH2*[fNcentBins];
	61	fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
	62	fHistDeltaPtEmbPtProbe = new TH2*[fNcentBins];
	63
	64	for (Int_t i = 0; i < fNcentBins; i++) {
65	fHistEmbJetsPtArea[i] = 0;
66	fHistEmbJetsCorrPtArea[i] = 0;
67	fHistEmbPartPtvsJetPt[i] = 0;
68	fHistEmbPartPtvsJetCorrPt[i] = 0;
69	fHistJetPtvsJetCorrPt[i] = 0;
70	fHistDistLeadPart2JetAxis[i] = 0;
71	fHistEmbBkgArea[i] = 0;
72	fHistRhoVSEmbBkg[i] = 0;
73	fHistDeltaPtEmbArea[i] = 0;
74	fHistDeltaPtEmbvsEP[i] = 0;
75	fHistDeltaPtEmbPtProbe[i] = 0;
76	}
77
78	SetMakeGeneralHistograms(kTRUE);
79	}
80
81	//________________________________________________________________________
82	AliAnalysisTaskDeltaPtJEmb::AliAnalysisTaskDeltaPtJEmb(const char *name) :
83	AliAnalysisTaskEmcalJet(name, kTRUE),
84	fJetsCont(0),
85	fTracksCont(0),
86	fCaloClustersCont(0),
87	fMinFractionShared(0.5),
88	fHistEmbJetsPtArea(0),
89	fHistEmbJetsCorrPtArea(0),
90	fHistEmbPartPtvsJetPt(0),
91	fHistEmbPartPtvsJetCorrPt(0),
92	fHistJetPtvsJetCorrPt(0),
93	fHistDistLeadPart2JetAxis(0),
94	fHistEmbBkgArea(0),
95	fHistRhoVSEmbBkg(0),
96	fHistDeltaPtEmbArea(0),
97	fHistDeltaPtEmbvsEP(0),
98	fHistDeltaPtEmbPtProbe(0),
99	fHistEmbJetsPhiEta(0),
100	fHistLeadPartPhiEta(0)
101	{
102	// Standard constructor.
103
104	fHistEmbJetsPtArea = new TH3*[fNcentBins];
105	fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
106	fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
107	fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
108	fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
109	fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
110	fHistEmbBkgArea = new TH2*[fNcentBins];
111	fHistRhoVSEmbBkg = new TH2*[fNcentBins];
112	fHistDeltaPtEmbArea = new TH2*[fNcentBins];
113	fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
114	fHistDeltaPtEmbPtProbe = new TH2*[fNcentBins];
115
116	for (Int_t i = 0; i < fNcentBins; i++) {
117	fHistEmbJetsPtArea[i] = 0;
118	fHistEmbJetsCorrPtArea[i] = 0;
119	fHistEmbPartPtvsJetPt[i] = 0;
120	fHistEmbPartPtvsJetCorrPt[i] = 0;
121	fHistJetPtvsJetCorrPt[i] = 0;
122	fHistDistLeadPart2JetAxis[i] = 0;
123	fHistEmbBkgArea[i] = 0;
124	fHistRhoVSEmbBkg[i] = 0;
125	fHistDeltaPtEmbArea[i] = 0;
126	fHistDeltaPtEmbvsEP[i] = 0;
127	fHistDeltaPtEmbPtProbe[i] = 0;
128	}
129
130	SetMakeGeneralHistograms(kTRUE);
131	}
132
133	//________________________________________________________________________
134	void AliAnalysisTaskDeltaPtJEmb::UserCreateOutputObjects()
135	{
136	// Create user output.
137
138	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
139
140	fJetsCont = GetJetContainer("Jets");
141	fTracksCont = GetParticleContainer("Tracks");
142	fCaloClustersCont = GetClusterContainer("CaloClusters");
143
144	fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
145	fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
146	fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
147	fOutput->Add(fHistEmbJetsPhiEta);
148
149	fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
150	fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
151	fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
152	fOutput->Add(fHistLeadPartPhiEta);
153
154	TString histname;
155
156	const Int_t nbinsZ = 12;
9dc6ae4b	157	Double_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
3a02e8fb	158
9dc6ae4b	159	Double_t *binsPt = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
	160	Double_t binsCorrPt = GenerateFixedBinArray(fNbins2, -fMaxBinPt, fMaxBinPt);
	161	Double_t *binsArea = GenerateFixedBinArray(50, 0, 2);
3a02e8fb	162
	163	for (Int_t i = 0; i < fNcentBins; i++) {
	164	histname = "fHistEmbJetsPtArea_";
	165	histname += i;
	166	fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
	167	fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
	168	fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
	169	fOutput->Add(fHistEmbJetsPtArea[i]);
	170
	171	histname = "fHistEmbJetsCorrPtArea_";
	172	histname += i;
	173	fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
	174	fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
	175	fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
	176	fOutput->Add(fHistEmbJetsCorrPtArea[i]);
	177
	178	histname = "fHistEmbPartPtvsJetPt_";
	179	histname += i;
	180	fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
	181	fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
	182	fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
	183	fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
	184	fOutput->Add(fHistEmbPartPtvsJetPt[i]);
	185
	186	histname = "fHistEmbPartPtvsJetCorrPt_";
	187	histname += i;
	188	fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
	189	fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
	190	fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
	191	fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
	192	fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
	193	fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
	194
	195	histname = "fHistJetPtvsJetCorrPt_";
	196	histname += i;
	197	fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
	198	fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
	199	fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
	200	fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
	201	fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
	202	fOutput->Add(fHistJetPtvsJetCorrPt[i]);
	203
	204	histname = "fHistDistLeadPart2JetAxis_";
	205	histname += i;
	206	fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
	207	fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
	208	fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
	209	fOutput->Add(fHistDistLeadPart2JetAxis[i]);
	210
	211	histname = "fHistEmbBkgArea_";
	212	histname += i;
	213	fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
	214	fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
	215	fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
	216	fOutput->Add(fHistEmbBkgArea[i]);
	217
	218	histname = "fHistRhoVSEmbBkg_";
	219	histname += i;
	220	fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
	221	fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
	222	fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
	223	fOutput->Add(fHistRhoVSEmbBkg[i]);
	224
	225	histname = "fHistDeltaPtEmbArea_";
226	histname += i;
227	fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(),
228	50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
229	fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
230	fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
231	fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
232	fOutput->Add(fHistDeltaPtEmbArea[i]);
233
234	histname = "fHistDeltaPtEmbvsEP_";
235	histname += i;
236	fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(),
237	402, -TMath::Pi()1.01, TMath::Pi()3.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
238	fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}");
239	fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
240	fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts");
241	fOutput->Add(fHistDeltaPtEmbvsEP[i]);
242
243	histname = "fHistDeltaPtEmbPtProbe_";
244	histname += i;
245	fHistDeltaPtEmbPtProbe[i] = new TH2F(histname.Data(), histname.Data(),
246	fNbins, fMinBinPt, fMaxBinPt, fNbins * 2, -fMaxBinPt, fMaxBinPt);
247	fHistDeltaPtEmbPtProbe[i]->GetXaxis()->SetTitle("#it{p}_{T,probe} (GeV/#it{c})");
248	fHistDeltaPtEmbPtProbe[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
249	fHistDeltaPtEmbPtProbe[i]->GetZaxis()->SetTitle("counts");
250	fOutput->Add(fHistDeltaPtEmbPtProbe[i]);
251	}
252
253	delete[] binsPt;
254	delete[] binsCorrPt;
255	delete[] binsArea;
256
257	PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
258	}
259
260	//________________________________________________________________________
261	Bool_t AliAnalysisTaskDeltaPtJEmb::FillHistograms()
262	{
263	// Fill histograms.
264
265	// ************
266	// Embedding
267	// Jets should already be matched, GetClosestJet gives matched partner. Cut on shared fraction to be applied
268	// _________________________________
269
270	if (fJetsCont) {
271	AliEmcalJet *jet1 = NULL;
272	AliEmcalJet *jet2 = NULL;
273	fJetsCont->ResetCurrentID();
274	while((jet1 = fJetsCont->GetNextAcceptJet())) {
275
276	jet2 = jet1->ClosestJet();
277	if(!jet2) continue;
278
279	Double_t fraction = fJetsCont->GetFractionSharedPt(jet1);
280	if(fMinFractionShared>0. && fraction<fMinFractionShared) continue;
281
282	TLorentzVector mom;
283	fJetsCont->GetLeadingHadronMomentum(mom,jet1);
284
285	Double_t distLeading2Jet = TMath::Sqrt((jet1->Eta() - mom.Eta()) * (jet1->Eta() - mom.Eta()) + (jet1->Phi() - mom.Phi()+TMath::Pi()) * (jet1->Phi() - mom.Phi()+TMath::Pi()));
286
287	fHistEmbPartPtvsJetPt[fCentBin]->Fill(jet2->Pt(), jet1->Pt());
288	fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(jet2->Pt(), jet1->Pt() - jet1->Area() * fRhoVal);
289	fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi()+TMath::Pi());
290	fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
291
292	fHistEmbJetsPtArea[fCentBin]->Fill(jet1->Area(), jet1->Pt(), mom.Pt());
293	fHistEmbJetsCorrPtArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area(), mom.Pt());
294	fHistEmbJetsPhiEta->Fill(jet1->Eta(), jet1->Phi());
295	fHistJetPtvsJetCorrPt[fCentBin]->Fill(jet1->Pt(), jet1->Pt() - fRhoVal * jet1->Area());
296
297	fHistEmbBkgArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - jet2->Pt());
298	fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * jet1->Area(), jet1->Pt() - jet2->Pt());
299	fHistDeltaPtEmbArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
300	fHistDeltaPtEmbvsEP[fCentBin]->Fill(jet1->Phi() - fEPV0, jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
301	fHistDeltaPtEmbPtProbe[fCentBin]->Fill(jet2->Pt(), jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
302	}
303	}
304
305	return kTRUE;
306	}
307
308	//________________________________________________________________________
309	void AliAnalysisTaskDeltaPtJEmb::ExecOnce()
310	{
311	// Initialize the analysis.
312
313	AliAnalysisTaskEmcalJet::ExecOnce();
314
315	if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
316	if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
317	if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
318	}