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

//
// Jet mass analysis task for jets recoiling from high pT hadron
//
// Author: M.Verweij

#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <THnSparse.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TProfile.h>
#include <TChain.h>
#include <TSystem.h>
#include <TFile.h>
#include <TKey.h>

#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliEmcalParticle.h"
#include "AliAnalysisManager.h"
#include "AliJetContainer.h"
#include "AliParticleContainer.h"

#include "AliAODEvent.h"

#include "AliAnalysisTaskEmcalHJetMass.h"

ClassImp(AliAnalysisTaskEmcalHJetMass)

//________________________________________________________________________
AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalHJetMass", kTRUE),
  fContainerBase(0),
  fContainerUnsub(1),
  fMinFractionShared(0),
  fUseUnsubJet(0),
  fJetMassType(kRaw),
  fDPhiHJetMax(0.6),
  fh1PtHadron(0),
  fh3PtJet1VsMassVsHPtAllSel(0),
  fh3PtJet1VsMassVsHPtTagged(0),
  fh3PtJet1VsMassVsHPtTaggedMatch(0),
  fh3PtJet1VsRatVsHPtAllSel(0),
  fh3PtJet1VsRatVsHPtTagged(0),
  fh3PtJet1VsRatVsHPtTaggedMatch(0)
{
  // Default constructor.

  fh1PtHadron                          = new TH1F*[fNcentBins];
  fh3PtJet1VsMassVsHPtAllSel        = new TH3F*[fNcentBins];
  fh3PtJet1VsMassVsHPtTagged        = new TH3F*[fNcentBins];
  fh3PtJet1VsMassVsHPtTaggedMatch   = new TH3F*[fNcentBins];
  fh3PtJet1VsRatVsHPtAllSel         = new TH3F*[fNcentBins];
  fh3PtJet1VsRatVsHPtTagged         = new TH3F*[fNcentBins];
  fh3PtJet1VsRatVsHPtTaggedMatch    = new TH3F*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh1PtHadron[i]                          = 0;
    fh3PtJet1VsMassVsHPtAllSel[i]        = 0;
    fh3PtJet1VsMassVsHPtTagged[i]        = 0;
    fh3PtJet1VsMassVsHPtTaggedMatch[i]   = 0;
    fh3PtJet1VsRatVsHPtAllSel[i]         = 0;
    fh3PtJet1VsRatVsHPtTagged[i]         = 0;
    fh3PtJet1VsRatVsHPtTaggedMatch[i]    = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),  
  fContainerBase(0),
  fContainerUnsub(1),
  fMinFractionShared(0),
  fUseUnsubJet(0),
  fJetMassType(kRaw),
  fDPhiHJetMax(0.6),
  fh1PtHadron(0),
  fh3PtJet1VsMassVsHPtAllSel(0),
  fh3PtJet1VsMassVsHPtTagged(0),
  fh3PtJet1VsMassVsHPtTaggedMatch(0),
  fh3PtJet1VsRatVsHPtAllSel(0),
  fh3PtJet1VsRatVsHPtTagged(0),
  fh3PtJet1VsRatVsHPtTaggedMatch(0)
{
  // Standard constructor.

  fh1PtHadron                          = new TH1F*[fNcentBins];
  fh3PtJet1VsMassVsHPtAllSel        = new TH3F*[fNcentBins];
  fh3PtJet1VsMassVsHPtTagged        = new TH3F*[fNcentBins];
  fh3PtJet1VsMassVsHPtTaggedMatch   = new TH3F*[fNcentBins];
  fh3PtJet1VsRatVsHPtAllSel         = new TH3F*[fNcentBins];
  fh3PtJet1VsRatVsHPtTagged         = new TH3F*[fNcentBins];
  fh3PtJet1VsRatVsHPtTaggedMatch    = new TH3F*[fNcentBins];
 
  for (Int_t i = 0; i < fNcentBins; i++) {
    fh1PtHadron[i]                          = 0;
    fh3PtJet1VsMassVsHPtAllSel[i]        = 0;
    fh3PtJet1VsMassVsHPtTagged[i]        = 0;
    fh3PtJet1VsMassVsHPtTaggedMatch[i]   = 0;
    fh3PtJet1VsRatVsHPtAllSel[i]         = 0;
    fh3PtJet1VsRatVsHPtTagged[i]         = 0;
    fh3PtJet1VsRatVsHPtTaggedMatch[i]    = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

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

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

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  const Int_t nBinsPt  = 200;
  const Double_t minPt = -50.;
  const Double_t maxPt = 150.;

  const Int_t nBinsM  = 100;
  const Double_t minM = -20.;
  const Double_t maxM = 80.;

  const Int_t nBinsR  = 100;
  const Double_t minR = -0.2;
  const Double_t maxR = 0.8;

  const Int_t nBinsPtH  = 100;
  const Double_t minPtH = 0.;
  const Double_t maxPtH = 100.;

  TString histName = "";
  TString histTitle = "";
  for (Int_t i = 0; i < fNcentBins; i++) {
    histName = TString::Format("fh1PtHadron_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,h}",histName.Data());
    fh1PtHadron[i] = new TH1F(histName.Data(),histTitle.Data(),200.,0.,200.);
    fOutput->Add(fh1PtHadron[i]);

    histName = TString::Format("fh3PtJet1VsMassVsHPtAllSel_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
    fh3PtJet1VsMassVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
    fOutput->Add(fh3PtJet1VsMassVsHPtAllSel[i]);

    histName = TString::Format("fh3PtJet1VsMassVsHPtTagged_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
    fh3PtJet1VsMassVsHPtTagged[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
    fOutput->Add(fh3PtJet1VsMassVsHPtTagged[i]);

    histName = TString::Format("fh3PtJet1VsMassVsHPtTaggedMatch_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
    fh3PtJet1VsMassVsHPtTaggedMatch[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
    fOutput->Add(fh3PtJet1VsMassVsHPtTaggedMatch[i]);

    //
    histName = TString::Format("fh3PtJet1VsRatVsHPtAllSel_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
    fh3PtJet1VsRatVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
    fOutput->Add(fh3PtJet1VsRatVsHPtAllSel[i]);

    histName = TString::Format("fh3PtJet1VsRatVsHPtTagged_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
    fh3PtJet1VsRatVsHPtTagged[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
    fOutput->Add(fh3PtJet1VsRatVsHPtTagged[i]);

    histName = TString::Format("fh3PtJet1VsRatVsHPtTaggedMatch_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
    fh3PtJet1VsRatVsHPtTaggedMatch[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
    fOutput->Add(fh3PtJet1VsRatVsHPtTaggedMatch[i]);
  }

  TH1::AddDirectory(oldStatus);

  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalHJetMass::Run()
{
  // Run analysis code here, if needed. It will be executed before FillHistograms().

  AliVParticle *vp = NULL;
  AliParticleContainer *pCont = GetParticleContainer(0);
  AliJetContainer      *jCont = GetJetContainer(fContainerBase);
  if(pCont) {
    pCont->ResetCurrentID();
    while((vp = pCont->GetNextAcceptParticle())) {
      fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all hadrons
      AliEmcalJet* jet = NULL;
      if(jCont) {
	jCont->ResetCurrentID();
	while((jet = jCont->GetNextAcceptJet())) {
	  Double_t dphi = GetDeltaPhi(vp,jet)+TMath::Pi();
	  if(dphi>fDPhiHJetMax) continue;
	  FillHJetHistograms(vp->Pt(),jet);
	}
      }
    }
  }
  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalHJetMass::FillHJetHistograms(Double_t pt, const AliEmcalJet *jet)
{
  // Fill hadron-jet histograms.
  Double_t ptJet = jet->Pt() - GetRhoVal(fContainerBase)*jet->Area();
  Double_t mJet = GetJetMass(jet);
  Double_t rat = -1.;
  if(ptJet<0. || ptJet>0.) rat = mJet/ptJet;

  fh3PtJet1VsMassVsHPtAllSel[fCentBin]->Fill(ptJet,mJet,pt);
  fh3PtJet1VsRatVsHPtAllSel[fCentBin]->Fill(ptJet,rat,pt);

  if(jet->GetTagStatus()<1 || !jet->GetTaggedJet())
    return kFALSE;
  fh3PtJet1VsMassVsHPtTagged[fCentBin]->Fill(ptJet,mJet,pt);
  fh3PtJet1VsRatVsHPtTagged[fCentBin]->Fill(ptJet,rat,pt);

  Double_t fraction = -1.;
  if(fUseUnsubJet) {
    AliEmcalJet *jetUS = NULL;
    AliJetContainer *jetContUS = GetJetContainer(fContainerUnsub);
    Int_t ifound = 0;
    Int_t ilab = -1;
    for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
      jetUS = jetContUS->GetJet(i);
      if(jetUS->GetLabel()==jet->GetLabel()) {
	ifound++;
	if(ifound==1) ilab = i;
      }
    }
    if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
    if(ifound==0) jetUS = 0x0;
    else          jetUS = jetContUS->GetJet(ilab);
    fraction = jetContUS->GetFractionSharedPt(jetUS);
  } else {
    AliJetContainer *jetCont = GetJetContainer(fContainerBase);
    fraction = jetCont->GetFractionSharedPt(jet);
  }
  
  if(fMinFractionShared>0. && fraction>fMinFractionShared) {
    fh3PtJet1VsMassVsHPtTaggedMatch[fCentBin]->Fill(ptJet,mJet,pt);
    fh3PtJet1VsRatVsHPtTaggedMatch[fCentBin]->Fill(ptJet,rat,pt);
  }
  return kTRUE;
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalHJetMass::GetJetMass(const AliEmcalJet *jet) const {
  //calc subtracted jet mass
  if(fJetMassType==kRaw)
    return jet->M();
  else if(fJetMassType==kDeriv)
    return jet->GetSecondOrderSubtracted();
  
  return 0;
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(const AliVParticle *vp, const AliEmcalJet* jet) const {
  // Calculate azimuthal angle between particle and jet. range:[-0.5\pi,1.5\pi]
  return GetDeltaPhi(vp->Phi(),jet->Phi());
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(Double_t phi1, Double_t phi2) const {
  // Calculate azimuthal angle between phi1 and phi2. range:[-0.5\pi,1.5\pi]
  Double_t dPhi = phi1-phi2;
  if(dPhi <-0.5*TMath::Pi())  dPhi += TMath::TwoPi();
  if(dPhi > 1.5*TMath::Pi())  dPhi -= TMath::TwoPi();
  return dPhi;
}


//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalHJetMass::RetrieveEventObjects() {
  //
  // retrieve event objects
  //
  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
    return kFALSE;

  return kTRUE;
}

//_______________________________________________________________________
void AliAnalysisTaskEmcalHJetMass::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
1705bab1	1	//
	2	// Jet mass analysis task for jets recoiling from high pT hadron
	3	//
	4	// Author: M.Verweij
	5
	6	#include <TClonesArray.h>
	7	#include <TH1F.h>
	8	#include <TH2F.h>
	9	#include <TH3F.h>
	10	#include <THnSparse.h>
	11	#include <TList.h>
	12	#include <TLorentzVector.h>
	13	#include <TProfile.h>
	14	#include <TChain.h>
	15	#include <TSystem.h>
	16	#include <TFile.h>
	17	#include <TKey.h>
	18
	19	#include "AliVCluster.h"
	20	#include "AliVTrack.h"
	21	#include "AliEmcalJet.h"
	22	#include "AliRhoParameter.h"
	23	#include "AliLog.h"
	24	#include "AliEmcalParticle.h"
	25	#include "AliAnalysisManager.h"
	26	#include "AliJetContainer.h"
	27	#include "AliParticleContainer.h"
	28
	29	#include "AliAODEvent.h"
	30
	31	#include "AliAnalysisTaskEmcalHJetMass.h"
	32
	33	ClassImp(AliAnalysisTaskEmcalHJetMass)
	34
	35	//________________________________________________________________________
	36	AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass() :
	37	AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalHJetMass", kTRUE),
	38	fContainerBase(0),
	39	fContainerUnsub(1),
	40	fMinFractionShared(0),
	41	fUseUnsubJet(0),
	42	fJetMassType(kRaw),
	43	fDPhiHJetMax(0.6),
	44	fh1PtHadron(0),
	45	fh3PtJet1VsMassVsHPtAllSel(0),
	46	fh3PtJet1VsMassVsHPtTagged(0),
	47	fh3PtJet1VsMassVsHPtTaggedMatch(0),
	48	fh3PtJet1VsRatVsHPtAllSel(0),
	49	fh3PtJet1VsRatVsHPtTagged(0),
	50	fh3PtJet1VsRatVsHPtTaggedMatch(0)
	51	{
	52	// Default constructor.
	53
	54	fh1PtHadron = new TH1F*[fNcentBins];
	55	fh3PtJet1VsMassVsHPtAllSel = new TH3F*[fNcentBins];
	56	fh3PtJet1VsMassVsHPtTagged = new TH3F*[fNcentBins];
	57	fh3PtJet1VsMassVsHPtTaggedMatch = new TH3F*[fNcentBins];
	58	fh3PtJet1VsRatVsHPtAllSel = new TH3F*[fNcentBins];
	59	fh3PtJet1VsRatVsHPtTagged = new TH3F*[fNcentBins];
	60	fh3PtJet1VsRatVsHPtTaggedMatch = new TH3F*[fNcentBins];
	61
	62	for (Int_t i = 0; i < fNcentBins; i++) {
	63	fh1PtHadron[i] = 0;
	64	fh3PtJet1VsMassVsHPtAllSel[i] = 0;
65	fh3PtJet1VsMassVsHPtTagged[i] = 0;
66	fh3PtJet1VsMassVsHPtTaggedMatch[i] = 0;
67	fh3PtJet1VsRatVsHPtAllSel[i] = 0;
68	fh3PtJet1VsRatVsHPtTagged[i] = 0;
69	fh3PtJet1VsRatVsHPtTaggedMatch[i] = 0;
70	}
71
72	SetMakeGeneralHistograms(kTRUE);
73	}
74
75	//________________________________________________________________________
76	AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass(const char *name) :
77	AliAnalysisTaskEmcalJet(name, kTRUE),
78	fContainerBase(0),
79	fContainerUnsub(1),
80	fMinFractionShared(0),
81	fUseUnsubJet(0),
82	fJetMassType(kRaw),
83	fDPhiHJetMax(0.6),
84	fh1PtHadron(0),
85	fh3PtJet1VsMassVsHPtAllSel(0),
86	fh3PtJet1VsMassVsHPtTagged(0),
87	fh3PtJet1VsMassVsHPtTaggedMatch(0),
88	fh3PtJet1VsRatVsHPtAllSel(0),
89	fh3PtJet1VsRatVsHPtTagged(0),
90	fh3PtJet1VsRatVsHPtTaggedMatch(0)
91	{
92	// Standard constructor.
93
94	fh1PtHadron = new TH1F*[fNcentBins];
95	fh3PtJet1VsMassVsHPtAllSel = new TH3F*[fNcentBins];
96	fh3PtJet1VsMassVsHPtTagged = new TH3F*[fNcentBins];
97	fh3PtJet1VsMassVsHPtTaggedMatch = new TH3F*[fNcentBins];
98	fh3PtJet1VsRatVsHPtAllSel = new TH3F*[fNcentBins];
99	fh3PtJet1VsRatVsHPtTagged = new TH3F*[fNcentBins];
100	fh3PtJet1VsRatVsHPtTaggedMatch = new TH3F*[fNcentBins];
101
102	for (Int_t i = 0; i < fNcentBins; i++) {
103	fh1PtHadron[i] = 0;
104	fh3PtJet1VsMassVsHPtAllSel[i] = 0;
105	fh3PtJet1VsMassVsHPtTagged[i] = 0;
106	fh3PtJet1VsMassVsHPtTaggedMatch[i] = 0;
107	fh3PtJet1VsRatVsHPtAllSel[i] = 0;
108	fh3PtJet1VsRatVsHPtTagged[i] = 0;
109	fh3PtJet1VsRatVsHPtTaggedMatch[i] = 0;
110	}
111
112	SetMakeGeneralHistograms(kTRUE);
113	}
114
115	//________________________________________________________________________
116	AliAnalysisTaskEmcalHJetMass::~AliAnalysisTaskEmcalHJetMass()
117	{
118	// Destructor.
119	}
120
121	//________________________________________________________________________
122	void AliAnalysisTaskEmcalHJetMass::UserCreateOutputObjects()
123	{
124	// Create user output.
125
126	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
127
128	Bool_t oldStatus = TH1::AddDirectoryStatus();
129	TH1::AddDirectory(kFALSE);
130
131	const Int_t nBinsPt = 200;
132	const Double_t minPt = -50.;
133	const Double_t maxPt = 150.;
134
135	const Int_t nBinsM = 100;
136	const Double_t minM = -20.;
137	const Double_t maxM = 80.;
138
139	const Int_t nBinsR = 100;
140	const Double_t minR = -0.2;
141	const Double_t maxR = 0.8;
142
143	const Int_t nBinsPtH = 100;
144	const Double_t minPtH = 0.;
145	const Double_t maxPtH = 100.;
146
147	TString histName = "";
148	TString histTitle = "";
149	for (Int_t i = 0; i < fNcentBins; i++) {
150	histName = TString::Format("fh1PtHadron_%d",i);
151	histTitle = TString::Format("%s;#it{p}_{T,h}",histName.Data());
152	fh1PtHadron[i] = new TH1F(histName.Data(),histTitle.Data(),200.,0.,200.);
153	fOutput->Add(fh1PtHadron[i]);
154
155	histName = TString::Format("fh3PtJet1VsMassVsHPtAllSel_%d",i);
156	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
157	fh3PtJet1VsMassVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
158	fOutput->Add(fh3PtJet1VsMassVsHPtAllSel[i]);
159
160	histName = TString::Format("fh3PtJet1VsMassVsHPtTagged_%d",i);
161	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
162	fh3PtJet1VsMassVsHPtTagged[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
163	fOutput->Add(fh3PtJet1VsMassVsHPtTagged[i]);
164
165	histName = TString::Format("fh3PtJet1VsMassVsHPtTaggedMatch_%d",i);
166	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
167	fh3PtJet1VsMassVsHPtTaggedMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
168	fOutput->Add(fh3PtJet1VsMassVsHPtTaggedMatch[i]);
169
170	//
171	histName = TString::Format("fh3PtJet1VsRatVsHPtAllSel_%d",i);
172	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
173	fh3PtJet1VsRatVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
174	fOutput->Add(fh3PtJet1VsRatVsHPtAllSel[i]);
175
176	histName = TString::Format("fh3PtJet1VsRatVsHPtTagged_%d",i);
177	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
178	fh3PtJet1VsRatVsHPtTagged[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
179	fOutput->Add(fh3PtJet1VsRatVsHPtTagged[i]);
180
181	histName = TString::Format("fh3PtJet1VsRatVsHPtTaggedMatch_%d",i);
182	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
183	fh3PtJet1VsRatVsHPtTaggedMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
184	fOutput->Add(fh3PtJet1VsRatVsHPtTaggedMatch[i]);
185	}
186
187	TH1::AddDirectory(oldStatus);
188
189	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
190	}
191
192	//________________________________________________________________________
193	Bool_t AliAnalysisTaskEmcalHJetMass::Run()
194	{
195	// Run analysis code here, if needed. It will be executed before FillHistograms().
196
197	AliVParticle *vp = NULL;
198	AliParticleContainer *pCont = GetParticleContainer(0);
199	AliJetContainer *jCont = GetJetContainer(fContainerBase);
200	if(pCont) {
201	pCont->ResetCurrentID();
202	while((vp = pCont->GetNextAcceptParticle())) {
203	fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all hadrons
204	AliEmcalJet* jet = NULL;
205	if(jCont) {
206	jCont->ResetCurrentID();
207	while((jet = jCont->GetNextAcceptJet())) {
208	Double_t dphi = GetDeltaPhi(vp,jet)+TMath::Pi();
209	if(dphi>fDPhiHJetMax) continue;
210	FillHJetHistograms(vp->Pt(),jet);
211	}
212	}
213	}
214	}
215	return kTRUE;
216	}
217
218	//________________________________________________________________________
219	Bool_t AliAnalysisTaskEmcalHJetMass::FillHJetHistograms(Double_t pt, const AliEmcalJet *jet)
220	{
221	// Fill hadron-jet histograms.
222	Double_t ptJet = jet->Pt() - GetRhoVal(fContainerBase)*jet->Area();
223	Double_t mJet = GetJetMass(jet);
224	Double_t rat = -1.;
225	if(ptJet<0. \|\| ptJet>0.) rat = mJet/ptJet;
226
227	fh3PtJet1VsMassVsHPtAllSel[fCentBin]->Fill(ptJet,mJet,pt);
228	fh3PtJet1VsRatVsHPtAllSel[fCentBin]->Fill(ptJet,rat,pt);
229
230	if(jet->GetTagStatus()<1 \|\| !jet->GetTaggedJet())
231	return kFALSE;
232	fh3PtJet1VsMassVsHPtTagged[fCentBin]->Fill(ptJet,mJet,pt);
233	fh3PtJet1VsRatVsHPtTagged[fCentBin]->Fill(ptJet,rat,pt);
234
235	Double_t fraction = -1.;
236	if(fUseUnsubJet) {
237	AliEmcalJet *jetUS = NULL;
238	AliJetContainer *jetContUS = GetJetContainer(fContainerUnsub);
239	Int_t ifound = 0;
240	Int_t ilab = -1;
241	for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
242	jetUS = jetContUS->GetJet(i);
243	if(jetUS->GetLabel()==jet->GetLabel()) {
244	ifound++;
245	if(ifound==1) ilab = i;
246	}
247	}
248	if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
249	if(ifound==0) jetUS = 0x0;
250	else jetUS = jetContUS->GetJet(ilab);
251	fraction = jetContUS->GetFractionSharedPt(jetUS);
252	} else {
253	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
254	fraction = jetCont->GetFractionSharedPt(jet);
255	}
256
257	if(fMinFractionShared>0. && fraction>fMinFractionShared) {
258	fh3PtJet1VsMassVsHPtTaggedMatch[fCentBin]->Fill(ptJet,mJet,pt);
259	fh3PtJet1VsRatVsHPtTaggedMatch[fCentBin]->Fill(ptJet,rat,pt);
260	}
261	return kTRUE;
262	}
263
264	//________________________________________________________________________
265	Double_t AliAnalysisTaskEmcalHJetMass::GetJetMass(const AliEmcalJet *jet) const {
266	//calc subtracted jet mass
267	if(fJetMassType==kRaw)
268	return jet->M();
269	else if(fJetMassType==kDeriv)
270	return jet->GetSecondOrderSubtracted();
271
272	return 0;
273	}
274
275	//________________________________________________________________________
276	Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(const AliVParticle vp, const AliEmcalJet jet) const {
277	// Calculate azimuthal angle between particle and jet. range:[-0.5\pi,1.5\pi]
278	return GetDeltaPhi(vp->Phi(),jet->Phi());
279	}
280
281	//________________________________________________________________________
282	Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(Double_t phi1, Double_t phi2) const {
283	// Calculate azimuthal angle between phi1 and phi2. range:[-0.5\pi,1.5\pi]
284	Double_t dPhi = phi1-phi2;
285	if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
286	if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
287	return dPhi;
288	}
289
290
291	//________________________________________________________________________
292	Bool_t AliAnalysisTaskEmcalHJetMass::RetrieveEventObjects() {
293	//
294	// retrieve event objects
295	//
296	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
297	return kFALSE;
298
299	return kTRUE;
300	}
301
302	//_______________________________________________________________________
303	void AliAnalysisTaskEmcalHJetMass::Terminate(Option_t *)
304	{
305	// Called once at the end of the analysis.
306	}
307