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

//
// Jet mass analysis task.
//
// 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 "AliMCEvent.h"
#include "AliGenPythiaEventHeader.h"
#include "AliAODMCHeader.h"
#include "AliMCEvent.h"
#include "AliAnalysisManager.h"
#include "AliJetContainer.h"

#include "AliAODEvent.h"

#include "AliAnalysisTaskEmcalJetMass.h"

ClassImp(AliAnalysisTaskEmcalJetMass)

//________________________________________________________________________
AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMass", kTRUE),
  fContainerBase(0),
  fMinFractionShared(0),
  fJetMassType(kRaw),
  fh3PtJet1VsMassVsLeadPtAllSel(0),
  fh3PtJet1VsMassVsLeadPtTagged(0),
  fh3PtJet1VsMassVsLeadPtTaggedMatch(0),
  fpPtVsMassJet1All(0),
  fpPtVsMassJet1Tagged(0),
  fpPtVsMassJet1TaggedMatch(0),
  fh2MassVsAreaJet1All(0),
  fh2MassVsAreaJet1Tagged(0),
  fh2MassVsAreaJet1TaggedMatch(0),
  fh2MassVsNConstJet1All(0),
  fh2MassVsNConstJet1Tagged(0),
  fh2MassVsNConstJet1TaggedMatch(0),
  fh2EtMassOverEtRSq(0)
{
  // Default constructor.

  fh3PtJet1VsMassVsLeadPtAllSel        = new TH3F*[fNcentBins];
  fh3PtJet1VsMassVsLeadPtTagged        = new TH3F*[fNcentBins];
  fh3PtJet1VsMassVsLeadPtTaggedMatch   = new TH3F*[fNcentBins];
  fpPtVsMassJet1All                    = new TProfile*[fNcentBins];
  fpPtVsMassJet1Tagged                 = new TProfile*[fNcentBins];
  fpPtVsMassJet1TaggedMatch            = new TProfile*[fNcentBins];
  fh2MassVsAreaJet1All                 = new TH2F*[fNcentBins];
  fh2MassVsAreaJet1Tagged              = new TH2F*[fNcentBins];
  fh2MassVsAreaJet1TaggedMatch         = new TH2F*[fNcentBins];
  fh2MassVsNConstJet1All               = new TH2F*[fNcentBins];
  fh2MassVsNConstJet1Tagged            = new TH2F*[fNcentBins];
  fh2MassVsNConstJet1TaggedMatch       = new TH2F*[fNcentBins];
  fh2EtMassOverEtRSq                   = new TH2F*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh3PtJet1VsMassVsLeadPtAllSel[i]        = 0;
    fh3PtJet1VsMassVsLeadPtTagged[i]        = 0;
    fh3PtJet1VsMassVsLeadPtTaggedMatch[i]   = 0;
    fpPtVsMassJet1All[i]                    = 0;
    fpPtVsMassJet1Tagged[i]                 = 0;
    fpPtVsMassJet1TaggedMatch[i]            = 0;
    fh2MassVsAreaJet1All[i]                 = 0;
    fh2MassVsAreaJet1Tagged[i]              = 0;
    fh2MassVsAreaJet1TaggedMatch[i]         = 0;
    fh2MassVsNConstJet1All[i]               = 0;
    fh2MassVsNConstJet1Tagged[i]            = 0;
    fh2MassVsNConstJet1TaggedMatch[i]       = 0;
    fh2EtMassOverEtRSq[i]                   = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),  
  fContainerBase(0),
  fMinFractionShared(0),
  fJetMassType(kRaw),
  fh3PtJet1VsMassVsLeadPtAllSel(0),
  fh3PtJet1VsMassVsLeadPtTagged(0),
  fh3PtJet1VsMassVsLeadPtTaggedMatch(0),
  fpPtVsMassJet1All(0),
  fpPtVsMassJet1Tagged(0),
  fpPtVsMassJet1TaggedMatch(0),
  fh2MassVsAreaJet1All(0),
  fh2MassVsAreaJet1Tagged(0),
  fh2MassVsAreaJet1TaggedMatch(0),
  fh2MassVsNConstJet1All(0),
  fh2MassVsNConstJet1Tagged(0),
  fh2MassVsNConstJet1TaggedMatch(0),
  fh2EtMassOverEtRSq(0)
{
  // Standard constructor.

  fh3PtJet1VsMassVsLeadPtAllSel        = new TH3F*[fNcentBins];
  fh3PtJet1VsMassVsLeadPtTagged        = new TH3F*[fNcentBins];
  fh3PtJet1VsMassVsLeadPtTaggedMatch   = new TH3F*[fNcentBins];
  fpPtVsMassJet1All                    = new TProfile*[fNcentBins];
  fpPtVsMassJet1Tagged                 = new TProfile*[fNcentBins];
  fpPtVsMassJet1TaggedMatch            = new TProfile*[fNcentBins];
  fh2MassVsAreaJet1All                 = new TH2F*[fNcentBins];
  fh2MassVsAreaJet1Tagged              = new TH2F*[fNcentBins];
  fh2MassVsAreaJet1TaggedMatch         = new TH2F*[fNcentBins];
  fh2MassVsNConstJet1All               = new TH2F*[fNcentBins];
  fh2MassVsNConstJet1Tagged            = new TH2F*[fNcentBins];
  fh2MassVsNConstJet1TaggedMatch       = new TH2F*[fNcentBins];
  fh2EtMassOverEtRSq                   = new TH2F*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh3PtJet1VsMassVsLeadPtAllSel[i]        = 0;
    fh3PtJet1VsMassVsLeadPtTagged[i]        = 0;
    fh3PtJet1VsMassVsLeadPtTaggedMatch[i]   = 0;
    fpPtVsMassJet1All[i]                    = 0;
    fpPtVsMassJet1Tagged[i]                 = 0;
    fpPtVsMassJet1TaggedMatch[i]            = 0;
    fh2MassVsAreaJet1All[i]                 = 0;
    fh2MassVsAreaJet1Tagged[i]              = 0;
    fh2MassVsAreaJet1TaggedMatch[i]         = 0;
    fh2MassVsNConstJet1All[i]               = 0;
    fh2MassVsNConstJet1Tagged[i]            = 0;
    fh2MassVsNConstJet1TaggedMatch[i]       = 0;
    fh2EtMassOverEtRSq[i]                   = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

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

//________________________________________________________________________
void AliAnalysisTaskEmcalJetMass::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 nBinsArea = 50;
  const Double_t minArea = 0.;
  const Double_t maxArea = 1.;

  const Int_t nBinsNConst = 100;
  const Double_t minNConst = 0.;
  const Double_t maxNConst = 500.;

  TString histName = "";
  TString histTitle = "";
  for (Int_t i = 0; i < fNcentBins; i++) {
    histName = TString::Format("fh3PtJet1VsMassVsLeadPtAllSel_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
    fh3PtJet1VsMassVsLeadPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
    fOutput->Add(fh3PtJet1VsMassVsLeadPtAllSel[i]);

    histName = TString::Format("fh3PtJet1VsMassVsLeadPtTagged_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
    fh3PtJet1VsMassVsLeadPtTagged[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
    fOutput->Add(fh3PtJet1VsMassVsLeadPtTagged[i]);

    histName = TString::Format("fh3PtJet1VsMassVsLeadPtTaggedMatch_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
    fh3PtJet1VsMassVsLeadPtTaggedMatch[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
    fOutput->Add(fh3PtJet1VsMassVsLeadPtTaggedMatch[i]);

    histName = TString::Format("fpPtVsMassJet1All_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
    fpPtVsMassJet1All[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
    fOutput->Add(fpPtVsMassJet1All[i]);

    histName = TString::Format("fpPtVsMassJet1Tagged_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
    fpPtVsMassJet1Tagged[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
    fOutput->Add(fpPtVsMassJet1Tagged[i]);

    histName = TString::Format("fpPtVsMassJet1TaggedMatch_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
    fpPtVsMassJet1TaggedMatch[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
    fOutput->Add(fpPtVsMassJet1TaggedMatch[i]);

    histName = TString::Format("fh2MassVsAreaJet1All_%d",i);
    histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
    fh2MassVsAreaJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
    fOutput->Add(fh2MassVsAreaJet1All[i]);

    histName = TString::Format("fh2MassVsAreaJet1Tagged_%d",i);
    histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
    fh2MassVsAreaJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
    fOutput->Add(fh2MassVsAreaJet1Tagged[i]);

    histName = TString::Format("fh2MassVsAreaJet1TaggedMatch_%d",i);
    histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
    fh2MassVsAreaJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
    fOutput->Add(fh2MassVsAreaJet1TaggedMatch[i]);

    histName = TString::Format("fh2MassVsNConstJet1All_%d",i);
    histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
    fh2MassVsNConstJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsNConst,minNConst,maxNConst);
    fOutput->Add(fh2MassVsNConstJet1All[i]);

    histName = TString::Format("fh2MassVsNConstJet1Tagged_%d",i);
    histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
    fh2MassVsNConstJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
    fOutput->Add(fh2MassVsNConstJet1Tagged[i]);

    histName = TString::Format("fh2MassVsNConstJet1TaggedMatch_%d",i);
    histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
    fh2MassVsNConstJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
    fOutput->Add(fh2MassVsNConstJet1TaggedMatch[i]);

    histName = TString::Format("fh2EtMassOverEtRSq_%d",i);
    histTitle = TString::Format("%s;#it{E}_{T};(#it{M}/(#it{E}_{T}#it{R}))^{2}",histName.Data());
    fh2EtMassOverEtRSq[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,1.);
    fOutput->Add(fh2EtMassOverEtRSq[i]);
  }

  // =========== Switch on Sumw2 for all histos ===========
  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
    TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
    if (h1){
      h1->Sumw2();
      continue;
    }
    THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
    if(hn)hn->Sumw2();
  }

  TH1::AddDirectory(oldStatus);

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

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

  return kTRUE;
}

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

  AliEmcalJet* jet1 = NULL;
  AliJetContainer *jetCont = GetJetContainer(fContainerBase);
  if(jetCont) {
    jetCont->ResetCurrentID();
    while((jet1 = jetCont->GetNextAcceptJet())) {

      Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
      Double_t mJet1 = GetJetMass(jet1);
      fh3PtJet1VsMassVsLeadPtAllSel[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
      fpPtVsMassJet1All[fCentBin]->Fill(ptJet1,mJet1);
      fh2MassVsAreaJet1All[fCentBin]->Fill(mJet1,jet1->Area());
      fh2MassVsNConstJet1All[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
      
      if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
	continue;

      fh3PtJet1VsMassVsLeadPtTagged[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
      fpPtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,mJet1);
      fh2MassVsAreaJet1Tagged[fCentBin]->Fill(mJet1,jet1->Area());
      fh2MassVsNConstJet1Tagged[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());

      Double_t fraction = jetCont->GetFractionSharedPt(jet1);
      if(fMinFractionShared>0. && fraction>fMinFractionShared) {
	fh3PtJet1VsMassVsLeadPtTaggedMatch[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
	fpPtVsMassJet1TaggedMatch[fCentBin]->Fill(ptJet1,mJet1);
	fh2MassVsAreaJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->Area());
	fh2MassVsNConstJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
      }
      
      Double_t Et2 = mJet1*mJet1 + ptJet1*ptJet1;
      Double_t Et = 0.;    Double_t massOverEtR = 0.;
      if(Et2>0.) Et = TMath::Sqrt(Et2);
      if((Et*jetCont->GetJetRadius())>0.) 
	massOverEtR = mJet1/(Et*jetCont->GetJetRadius());
      fh2EtMassOverEtRSq[fCentBin]->Fill(Et,massOverEtR*massOverEtR);
    }
  }
  
  return kTRUE;
}

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

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

  return kTRUE;
}

//_______________________________________________________________________
void AliAnalysisTaskEmcalJetMass::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
8612dfc8	1	//
	2	// Jet mass analysis task.
	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 "AliMCEvent.h"
	26	#include "AliGenPythiaEventHeader.h"
	27	#include "AliAODMCHeader.h"
	28	#include "AliMCEvent.h"
	29	#include "AliAnalysisManager.h"
	30	#include "AliJetContainer.h"
	31
	32	#include "AliAODEvent.h"
	33
	34	#include "AliAnalysisTaskEmcalJetMass.h"
	35
	36	ClassImp(AliAnalysisTaskEmcalJetMass)
	37
	38	//________________________________________________________________________
	39	AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass() :
	40	AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMass", kTRUE),
	41	fContainerBase(0),
ba7663ad	42	fMinFractionShared(0),
f831dc9f	43	fJetMassType(kRaw),
713c5683	44	fh3PtJet1VsMassVsLeadPtAllSel(0),
	45	fh3PtJet1VsMassVsLeadPtTagged(0),
	46	fh3PtJet1VsMassVsLeadPtTaggedMatch(0),
8612dfc8	47	fpPtVsMassJet1All(0),
8612dfc8	48	fpPtVsMassJet1Tagged(0),
4d3b366f	49	fpPtVsMassJet1TaggedMatch(0),
8612dfc8	50	fh2MassVsAreaJet1All(0),
02e33ee8	51	fh2MassVsAreaJet1Tagged(0),
4d3b366f	52	fh2MassVsAreaJet1TaggedMatch(0),
ba7663ad	53	fh2MassVsNConstJet1All(0),
ba7663ad	54	fh2MassVsNConstJet1Tagged(0),
4d3b366f	55	fh2MassVsNConstJet1TaggedMatch(0),
02e33ee8	56	fh2EtMassOverEtRSq(0)
8612dfc8	57	{
	58	// Default constructor.
	59
713c5683	60	fh3PtJet1VsMassVsLeadPtAllSel = new TH3F*[fNcentBins];
	61	fh3PtJet1VsMassVsLeadPtTagged = new TH3F*[fNcentBins];
	62	fh3PtJet1VsMassVsLeadPtTaggedMatch = new TH3F*[fNcentBins];
	63	fpPtVsMassJet1All = new TProfile*[fNcentBins];
	64	fpPtVsMassJet1Tagged = new TProfile*[fNcentBins];
	65	fpPtVsMassJet1TaggedMatch = new TProfile*[fNcentBins];
	66	fh2MassVsAreaJet1All = new TH2F*[fNcentBins];
	67	fh2MassVsAreaJet1Tagged = new TH2F*[fNcentBins];
	68	fh2MassVsAreaJet1TaggedMatch = new TH2F*[fNcentBins];
	69	fh2MassVsNConstJet1All = new TH2F*[fNcentBins];
	70	fh2MassVsNConstJet1Tagged = new TH2F*[fNcentBins];
	71	fh2MassVsNConstJet1TaggedMatch = new TH2F*[fNcentBins];
	72	fh2EtMassOverEtRSq = new TH2F*[fNcentBins];
8612dfc8	73
8612dfc8	74	for (Int_t i = 0; i < fNcentBins; i++) {
713c5683	75	fh3PtJet1VsMassVsLeadPtAllSel[i] = 0;
	76	fh3PtJet1VsMassVsLeadPtTagged[i] = 0;
	77	fh3PtJet1VsMassVsLeadPtTaggedMatch[i] = 0;
	78	fpPtVsMassJet1All[i] = 0;
	79	fpPtVsMassJet1Tagged[i] = 0;
	80	fpPtVsMassJet1TaggedMatch[i] = 0;
	81	fh2MassVsAreaJet1All[i] = 0;
	82	fh2MassVsAreaJet1Tagged[i] = 0;
	83	fh2MassVsAreaJet1TaggedMatch[i] = 0;
	84	fh2MassVsNConstJet1All[i] = 0;
	85	fh2MassVsNConstJet1Tagged[i] = 0;
	86	fh2MassVsNConstJet1TaggedMatch[i] = 0;
	87	fh2EtMassOverEtRSq[i] = 0;
8612dfc8	88	}
	89
	90	SetMakeGeneralHistograms(kTRUE);
8612dfc8	91	}
	92
	93	//________________________________________________________________________
	94	AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) :
	95	AliAnalysisTaskEmcalJet(name, kTRUE),
	96	fContainerBase(0),
ba7663ad	97	fMinFractionShared(0),
f831dc9f	98	fJetMassType(kRaw),
713c5683	99	fh3PtJet1VsMassVsLeadPtAllSel(0),
	100	fh3PtJet1VsMassVsLeadPtTagged(0),
	101	fh3PtJet1VsMassVsLeadPtTaggedMatch(0),
8612dfc8	102	fpPtVsMassJet1All(0),
8612dfc8	103	fpPtVsMassJet1Tagged(0),
4d3b366f	104	fpPtVsMassJet1TaggedMatch(0),
8612dfc8	105	fh2MassVsAreaJet1All(0),
02e33ee8	106	fh2MassVsAreaJet1Tagged(0),
4d3b366f	107	fh2MassVsAreaJet1TaggedMatch(0),
ba7663ad	108	fh2MassVsNConstJet1All(0),
ba7663ad	109	fh2MassVsNConstJet1Tagged(0),
4d3b366f	110	fh2MassVsNConstJet1TaggedMatch(0),
02e33ee8	111	fh2EtMassOverEtRSq(0)
8612dfc8	112	{
	113	// Standard constructor.
	114
713c5683	115	fh3PtJet1VsMassVsLeadPtAllSel = new TH3F*[fNcentBins];
	116	fh3PtJet1VsMassVsLeadPtTagged = new TH3F*[fNcentBins];
	117	fh3PtJet1VsMassVsLeadPtTaggedMatch = new TH3F*[fNcentBins];
	118	fpPtVsMassJet1All = new TProfile*[fNcentBins];
	119	fpPtVsMassJet1Tagged = new TProfile*[fNcentBins];
	120	fpPtVsMassJet1TaggedMatch = new TProfile*[fNcentBins];
	121	fh2MassVsAreaJet1All = new TH2F*[fNcentBins];
	122	fh2MassVsAreaJet1Tagged = new TH2F*[fNcentBins];
	123	fh2MassVsAreaJet1TaggedMatch = new TH2F*[fNcentBins];
	124	fh2MassVsNConstJet1All = new TH2F*[fNcentBins];
	125	fh2MassVsNConstJet1Tagged = new TH2F*[fNcentBins];
	126	fh2MassVsNConstJet1TaggedMatch = new TH2F*[fNcentBins];
	127	fh2EtMassOverEtRSq = new TH2F*[fNcentBins];
8612dfc8	128
8612dfc8	129	for (Int_t i = 0; i < fNcentBins; i++) {
713c5683	130	fh3PtJet1VsMassVsLeadPtAllSel[i] = 0;
	131	fh3PtJet1VsMassVsLeadPtTagged[i] = 0;
	132	fh3PtJet1VsMassVsLeadPtTaggedMatch[i] = 0;
	133	fpPtVsMassJet1All[i] = 0;
	134	fpPtVsMassJet1Tagged[i] = 0;
	135	fpPtVsMassJet1TaggedMatch[i] = 0;
	136	fh2MassVsAreaJet1All[i] = 0;
	137	fh2MassVsAreaJet1Tagged[i] = 0;
	138	fh2MassVsAreaJet1TaggedMatch[i] = 0;
	139	fh2MassVsNConstJet1All[i] = 0;
	140	fh2MassVsNConstJet1Tagged[i] = 0;
	141	fh2MassVsNConstJet1TaggedMatch[i] = 0;
	142	fh2EtMassOverEtRSq[i] = 0;
8612dfc8	143	}
	144
	145	SetMakeGeneralHistograms(kTRUE);
	146	}
	147
	148	//________________________________________________________________________
	149	AliAnalysisTaskEmcalJetMass::~AliAnalysisTaskEmcalJetMass()
	150	{
	151	// Destructor.
	152	}
	153
	154	//________________________________________________________________________
	155	void AliAnalysisTaskEmcalJetMass::UserCreateOutputObjects()
	156	{
	157	// Create user output.
	158
	159	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
	160
	161	Bool_t oldStatus = TH1::AddDirectoryStatus();
	162	TH1::AddDirectory(kFALSE);
	163
713c5683	164	const Int_t nBinsPt = 200;
8612dfc8	165	const Double_t minPt = -50.;
713c5683	166	const Double_t maxPt = 150.;
8612dfc8	167
713c5683	168	const Int_t nBinsM = 100;
	169	const Double_t minM = -20.;
	170	const Double_t maxM = 80.;
372d124e	171
713c5683	172	const Int_t nBinsArea = 50;
8612dfc8	173	const Double_t minArea = 0.;
	174	const Double_t maxArea = 1.;
	175
ba7663ad	176	const Int_t nBinsNConst = 100;
	177	const Double_t minNConst = 0.;
	178	const Double_t maxNConst = 500.;
	179
8612dfc8	180	TString histName = "";
	181	TString histTitle = "";
	182	for (Int_t i = 0; i < fNcentBins; i++) {
713c5683	183	histName = TString::Format("fh3PtJet1VsMassVsLeadPtAllSel_%d",i);
	184	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
	185	fh3PtJet1VsMassVsLeadPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
	186	fOutput->Add(fh3PtJet1VsMassVsLeadPtAllSel[i]);
	187
	188	histName = TString::Format("fh3PtJet1VsMassVsLeadPtTagged_%d",i);
	189	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
	190	fh3PtJet1VsMassVsLeadPtTagged[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
	191	fOutput->Add(fh3PtJet1VsMassVsLeadPtTagged[i]);
	192
	193	histName = TString::Format("fh3PtJet1VsMassVsLeadPtTaggedMatch_%d",i);
	194	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
	195	fh3PtJet1VsMassVsLeadPtTaggedMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
	196	fOutput->Add(fh3PtJet1VsMassVsLeadPtTaggedMatch[i]);
4d3b366f	197
8612dfc8	198	histName = TString::Format("fpPtVsMassJet1All_%d",i);
8612dfc8	199	histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
713c5683	200	fpPtVsMassJet1All[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
8612dfc8	201	fOutput->Add(fpPtVsMassJet1All[i]);
	202
	203	histName = TString::Format("fpPtVsMassJet1Tagged_%d",i);
	204	histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
713c5683	205	fpPtVsMassJet1Tagged[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
8612dfc8	206	fOutput->Add(fpPtVsMassJet1Tagged[i]);
8612dfc8	207
4d3b366f	208	histName = TString::Format("fpPtVsMassJet1TaggedMatch_%d",i);
4d3b366f	209	histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
713c5683	210	fpPtVsMassJet1TaggedMatch[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
4d3b366f	211	fOutput->Add(fpPtVsMassJet1TaggedMatch[i]);
4d3b366f	212
8612dfc8	213	histName = TString::Format("fh2MassVsAreaJet1All_%d",i);
8612dfc8	214	histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
372d124e	215	fh2MassVsAreaJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
8612dfc8	216	fOutput->Add(fh2MassVsAreaJet1All[i]);
	217
	218	histName = TString::Format("fh2MassVsAreaJet1Tagged_%d",i);
	219	histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
372d124e	220	fh2MassVsAreaJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
8612dfc8	221	fOutput->Add(fh2MassVsAreaJet1Tagged[i]);
8612dfc8	222
4d3b366f	223	histName = TString::Format("fh2MassVsAreaJet1TaggedMatch_%d",i);
	224	histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
	225	fh2MassVsAreaJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
	226	fOutput->Add(fh2MassVsAreaJet1TaggedMatch[i]);
	227
ba7663ad	228	histName = TString::Format("fh2MassVsNConstJet1All_%d",i);
ba7663ad	229	histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
372d124e	230	fh2MassVsNConstJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsNConst,minNConst,maxNConst);
ba7663ad	231	fOutput->Add(fh2MassVsNConstJet1All[i]);
	232
	233	histName = TString::Format("fh2MassVsNConstJet1Tagged_%d",i);
	234	histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
	235	fh2MassVsNConstJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
	236	fOutput->Add(fh2MassVsNConstJet1Tagged[i]);
	237
4d3b366f	238	histName = TString::Format("fh2MassVsNConstJet1TaggedMatch_%d",i);
	239	histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
	240	fh2MassVsNConstJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
	241	fOutput->Add(fh2MassVsNConstJet1TaggedMatch[i]);
	242
02e33ee8	243	histName = TString::Format("fh2EtMassOverEtRSq_%d",i);
	244	histTitle = TString::Format("%s;#it{E}_{T};(#it{M}/(#it{E}_{T}#it{R}))^{2}",histName.Data());
	245	fh2EtMassOverEtRSq[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,1.);
	246	fOutput->Add(fh2EtMassOverEtRSq[i]);
8612dfc8	247	}
	248
	249	// =========== Switch on Sumw2 for all histos ===========
	250	for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
	251	TH1 h1 = dynamic_cast<TH1>(fOutput->At(i));
	252	if (h1){
	253	h1->Sumw2();
	254	continue;
	255	}
	256	THnSparse hn = dynamic_cast<THnSparse>(fOutput->At(i));
	257	if(hn)hn->Sumw2();
	258	}
	259
	260	TH1::AddDirectory(oldStatus);
	261
	262	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
	263	}
	264
	265	//________________________________________________________________________
	266	Bool_t AliAnalysisTaskEmcalJetMass::Run()
	267	{
	268	// Run analysis code here, if needed. It will be executed before FillHistograms().
	269
	270	return kTRUE;
	271	}
	272
	273	//________________________________________________________________________
	274	Bool_t AliAnalysisTaskEmcalJetMass::FillHistograms()
	275	{
	276	// Fill histograms.
	277
ba7663ad	278	AliEmcalJet* jet1 = NULL;
02e33ee8	279	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
	280	if(jetCont) {
	281	jetCont->ResetCurrentID();
	282	while((jet1 = jetCont->GetNextAcceptJet())) {
ba7663ad	283
923da004	284	Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
4d3b366f	285	Double_t mJet1 = GetJetMass(jet1);
713c5683	286	fh3PtJet1VsMassVsLeadPtAllSel[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
4d3b366f	287	fpPtVsMassJet1All[fCentBin]->Fill(ptJet1,mJet1);
	288	fh2MassVsAreaJet1All[fCentBin]->Fill(mJet1,jet1->Area());
	289	fh2MassVsNConstJet1All[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
02e33ee8	290
ba7663ad	291	if(jet1->GetTagStatus()<1 \|\| !jet1->GetTaggedJet())
02e33ee8	292	continue;
ba7663ad	293
713c5683	294	fh3PtJet1VsMassVsLeadPtTagged[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
4d3b366f	295	fpPtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,mJet1);
	296	fh2MassVsAreaJet1Tagged[fCentBin]->Fill(mJet1,jet1->Area());
	297	fh2MassVsNConstJet1Tagged[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
	298
	299	Double_t fraction = jetCont->GetFractionSharedPt(jet1);
	300	if(fMinFractionShared>0. && fraction>fMinFractionShared) {
713c5683	301	fh3PtJet1VsMassVsLeadPtTaggedMatch[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
4d3b366f	302	fpPtVsMassJet1TaggedMatch[fCentBin]->Fill(ptJet1,mJet1);
	303	fh2MassVsAreaJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->Area());
	304	fh2MassVsNConstJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
	305	}
02e33ee8	306
713c5683	307	Double_t Et2 = mJet1mJet1 + ptJet1ptJet1;
ef7f972e	308	Double_t Et = 0.; Double_t massOverEtR = 0.;
	309	if(Et2>0.) Et = TMath::Sqrt(Et2);
	310	if((Et*jetCont->GetJetRadius())>0.)
713c5683	311	massOverEtR = mJet1/(Et*jetCont->GetJetRadius());
ef7f972e	312	fh2EtMassOverEtRSq[fCentBin]->Fill(Et,massOverEtR*massOverEtR);
02e33ee8	313	}
8612dfc8	314	}
923da004	315
8612dfc8	316	return kTRUE;
8612dfc8	317	}
8612dfc8	318
2e31a073	319	//________________________________________________________________________
2e31a073	320	Double_t AliAnalysisTaskEmcalJetMass::GetJetMass(AliEmcalJet *jet) {
4d3b366f	321	//calc subtracted jet mass
f831dc9f	322	if(fJetMassType==kRaw)
	323	return jet->M();
	324	else if(fJetMassType==kDeriv)
	325	return jet->GetSecondOrderSubtracted();
	326
	327	return 0;
2e31a073	328	}
2e31a073	329
8612dfc8	330	//________________________________________________________________________
	331	Bool_t AliAnalysisTaskEmcalJetMass::RetrieveEventObjects() {
	332	//
	333	// retrieve event objects
	334	//
8612dfc8	335	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
	336	return kFALSE;
	337
	338	return kTRUE;
8612dfc8	339	}
	340
	341	//_______________________________________________________________________
	342	void AliAnalysisTaskEmcalJetMass::Terminate(Option_t *)
	343	{
	344	// Called once at the end of the analysis.
	345	}
	346