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

//
// Do subtraction for jet shapes using constituents arXiv:1403.3108
//
// Author: M.Verweij

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

#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliEmcalParticle.h"
#include "AliMCEvent.h"
#include "AliAODEvent.h"
#include "AliGenPythiaEventHeader.h"
#include "AliAODMCHeader.h"
#include "AliAnalysisManager.h"
#include "AliJetContainer.h"
#include "AliParticleContainer.h"

#include "AliAnalysisTaskJetShapeConst.h"

ClassImp(AliAnalysisTaskJetShapeConst)

//________________________________________________________________________
AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeConst", kTRUE),
  fContainerBase(0),
  fContainerSub(1),
  fContainerNoEmb(2),
  fMinFractionShared(0),
  fSingleTrackEmb(kFALSE),
  fCreateTree(kFALSE),
  fJetMassVarType(kMass),
  fUseSumw2(0),
  fTreeJetBkg(),
  fJet1Vec(new TLorentzVector()),
  fJet2Vec(new TLorentzVector()),
  fJetSubVec(new TLorentzVector()),
  fArea(0),
  fAreaPhi(0),
  fAreaEta(0),
  fRho(0),
  fRhoM(0),
  fNConst(0),
  fMatch(0),
  fh2MSubMatch(0x0),
  fh2MSubPtRawAll(0x0),
  fh3MSubPtRawDRMatch(0x0),
  fh3MSubPtTrueLeadPt(0x0),
  fh3MTruePtTrueLeadPt(0x0),
  fh3PtTrueDeltaMLeadPt(0x0),
  fh3PtTrueDeltaMRelLeadPt(0x0),
  fhnMassResponse(0x0)
{
  // Default constructor.

  fh2MSubMatch             = new TH2F*[fNcentBins];
  fh2MSubPtRawAll          = new TH2F*[fNcentBins];
  fh3MSubPtRawDRMatch      = new TH3F*[fNcentBins];
  fh3MSubPtTrueLeadPt      = new TH3F*[fNcentBins];
  fh3MTruePtTrueLeadPt     = new TH3F*[fNcentBins];
  fh3PtTrueDeltaMLeadPt    = new TH3F*[fNcentBins];
  fh3PtTrueDeltaMRelLeadPt = new TH3F*[fNcentBins];
  fhnMassResponse          = new THnSparse*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh2MSubMatch[i]             = 0;
    fh2MSubPtRawAll[i]          = 0;
    fh3MSubPtRawDRMatch[i]      = 0;
    fh3MSubPtTrueLeadPt[i]      = 0;
    fh3MTruePtTrueLeadPt[i]     = 0;
    fh3PtTrueDeltaMLeadPt[i]    = 0;
    fh3PtTrueDeltaMRelLeadPt[i] = 0;
    fhnMassResponse[i]          = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
  if(fCreateTree) DefineOutput(2, TTree::Class());
}

//________________________________________________________________________
AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),  
  fContainerBase(0),
  fContainerSub(1),
  fContainerNoEmb(2),
  fMinFractionShared(0),
  fSingleTrackEmb(kFALSE),
  fCreateTree(kFALSE),
  fJetMassVarType(kMass),
  fUseSumw2(0),
  fTreeJetBkg(0),
  fJet1Vec(new TLorentzVector()),
  fJet2Vec(new TLorentzVector()),
  fJetSubVec(new TLorentzVector()),
  fArea(0),
  fAreaPhi(0),
  fAreaEta(0),
  fRho(0),
  fRhoM(0),
  fNConst(0),
  fMatch(0),
  fh2MSubMatch(0x0),
  fh2MSubPtRawAll(0x0),
  fh3MSubPtRawDRMatch(0x0),
  fh3MSubPtTrueLeadPt(0x0),
  fh3MTruePtTrueLeadPt(0x0),
  fh3PtTrueDeltaMLeadPt(0x0),
  fh3PtTrueDeltaMRelLeadPt(0x0),
  fhnMassResponse(0x0)
{
  // Standard constructor.

  fh2MSubMatch             = new TH2F*[fNcentBins];
  fh2MSubPtRawAll          = new TH2F*[fNcentBins];
  fh3MSubPtRawDRMatch      = new TH3F*[fNcentBins];
  fh3MSubPtTrueLeadPt      = new TH3F*[fNcentBins];
  fh3MTruePtTrueLeadPt     = new TH3F*[fNcentBins];
  fh3PtTrueDeltaMLeadPt    = new TH3F*[fNcentBins];
  fh3PtTrueDeltaMRelLeadPt = new TH3F*[fNcentBins];
  fhnMassResponse          = new THnSparse*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh2MSubMatch[i]             = 0;
    fh2MSubPtRawAll[i]          = 0;
    fh3MSubPtRawDRMatch[i]      = 0;
    fh3MSubPtTrueLeadPt[i]      = 0;
    fh3MTruePtTrueLeadPt[i]     = 0;
    fh3PtTrueDeltaMLeadPt[i]    = 0;
    fh3PtTrueDeltaMRelLeadPt[i] = 0;
    fhnMassResponse[i]          = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
  if(fCreateTree) DefineOutput(2, TTree::Class());
}

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

//________________________________________________________________________
void AliAnalysisTaskJetShapeConst::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.;

  Int_t nBinsM  = 100;
  Double_t minM = -20.;
  Double_t maxM = 80.;
  if(fJetMassVarType==kRatMPt) {
    nBinsM = 100;
    minM   = -0.2;
    maxM   = 0.8;
  }

  Int_t nBinsDM  = 100;
  Double_t minDM = -25.;
  Double_t maxDM = 25.;
  if(fJetMassVarType==kRatMPt) {
    nBinsDM = 100;
    minDM   = -0.5;
    maxDM   = 0.5;
  }

  const Int_t nBinsDRToLJ  = 20; //distance to leading jet in Pb-Pb only event
  const Double_t minDRToLJ = 0.;
  const Double_t maxDRToLJ = 1.;

  const Int_t nBinsPtLead = 20;
  const Double_t minPtLead = 0.;
  const Double_t maxPtLead = 20.;

  //Binning for THnSparse
  const Int_t nBinsSparse0 = 5;
  const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsPtLead};
  const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt, minPtLead};
  const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt, maxPtLead};

  TString histName = "";
  TString histTitle = "";
  TString varName = "#it{M}_{jet}";
  if(fJetMassVarType==kRatMPt) varName = "#it{M}_{jet}/#it{p}_{T,jet}";

  for (Int_t i = 0; i < fNcentBins; i++) {
    histName = Form("fh2MSubMatch_%d",i);
    histTitle = Form("fh2MSubMatch_%d;%s;match",i,varName.Data());
    fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
    fOutput->Add(fh2MSubMatch[i]);

    histName = Form("fh2MSubPtRawAll_%d",i);
    histTitle = Form("fh2MSubPtRawAll_%d;%s;#it{p}_{T}",i,varName.Data());
    fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
    fOutput->Add(fh2MSubPtRawAll[i]);

    histName = Form("fh3MSubPtRawDRMatch_%d",i);
    histTitle = Form("fh3MSubPtRawDRMatch_%d;%s;#it{p}_{T}",i,varName.Data());
    fh3MSubPtRawDRMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
    fOutput->Add(fh3MSubPtRawDRMatch[i]);

    histName = Form("fh3MSubPtTrueLeadPt_%d",i);
    histTitle = Form("fh3MSubPtTrueLeadPt_%d;%s;#it{p}_{T}",i,varName.Data());
    fh3MSubPtTrueLeadPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsPtLead,minPtLead,maxPtLead);
    fOutput->Add(fh3MSubPtTrueLeadPt[i]);

    histName = Form("fh3MTruePtTrueLeadPt_%d",i);
    histTitle = Form("fh3MTruePtTrueLeadPt_%d;%s;#it{p}_{T}",i,varName.Data());
    fh3MTruePtTrueLeadPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsPtLead,minPtLead,maxPtLead);
    fOutput->Add(fh3MTruePtTrueLeadPt[i]);

    histName = Form("fh3PtTrueDeltaMLeadPt_%d",i);
    histTitle = Form("fh3PtTrueDeltaMLeadPt_%d;#it{p}_{T,true};#Delta %s",i,varName.Data());
    fh3PtTrueDeltaMLeadPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDM,minDM,maxDM,nBinsPtLead,minPtLead,maxPtLead);
    fOutput->Add(fh3PtTrueDeltaMLeadPt[i]);

    histName = Form("fh3PtTrueDeltaMRelLeadPt_%d",i);
    histTitle = Form("fh3PtTrueDeltaMRelLeadPt_%d;#it{p}_{T,true};Rel #Delta %s",i,varName.Data());
    fh3PtTrueDeltaMRelLeadPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,400,-1.,3.,nBinsPtLead,minPtLead,maxPtLead);
    fOutput->Add(fh3PtTrueDeltaMRelLeadPt[i]);

    histName = Form("fhnMassResponse_%d",i);
    histTitle = Form("fhnMassResponse_%d;%s sub;%s true;#it{p}_{T,sub};#it{p}_{T,true};#it{p}_{T,lead trk}",i,varName.Data(),varName.Data());
    fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
    fOutput->Add(fhnMassResponse[i]);
  }

  if(fUseSumw2) {
    // =========== 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);

  // Create a tree.
  if(fCreateTree) {
    fTreeJetBkg = new TTree("fTreeJetSubConst", "fTreeJetSubConst");
    fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
    fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
    fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
    fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
    fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
    fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
    fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
    fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
    fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
  }
  
  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
  if(fCreateTree) PostData(2, fTreeJetBkg);
}

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

  return kTRUE;
}

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

  AliEmcalJet* jet1  = NULL; //AA jet
  AliEmcalJet *jet2  = NULL; //Embedded Pythia jet
  //  AliEmcalJet *jet1T = NULL; //tagged AA jet
  //  AliEmcalJet *jet2T = NULL; //tagged Pythia jet
  AliEmcalJet *jetS = NULL;  //subtracted jet
  AliJetContainer *jetCont = GetJetContainer(fContainerBase);
  AliJetContainer *jetContS = GetJetContainer(fContainerSub);

  //Get leading jet in Pb-Pb event without embedded objects
  AliJetContainer *jetContNoEmb = GetJetContainer(fContainerNoEmb);
  AliEmcalJet *jetL = NULL;
  if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");

  AliDebug(11,Form("NJets  Incl: %d  Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
  if(jetCont) {
    jetCont->ResetCurrentID();
    while((jet1 = jetCont->GetNextAcceptJet())) {
      jet2  = NULL;
      jetS  = NULL;
      
      //Get constituent subtacted version of jet
      Int_t ifound = 0;
      Int_t ilab = -1;
      for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
	//if(ifound==1) continue;
	jetS = jetContS->GetJet(i);
	if(jetS->GetLabel()==jet1->GetLabel()) { // && jetS->Pt()>0.) {
	  ifound++;
	  if(ifound==1) ilab = i;
	}
      }
      if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
      if(ifound==0) jetS = 0x0;
      else jetS = jetContS->GetJet(ilab);

      Double_t mjetS = jetS->M();
      Double_t ptjet1 = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
      Double_t var = mjetS;
      if(fJetMassVarType==kRatMPt) {
	if(ptjet1>0. || ptjet1<0.) var = mjetS/ptjet1;
	else var = -999.;
      }

      //Fill histograms for all AA jets
      fh2MSubPtRawAll[fCentBin]->Fill(var,ptjet1);

      Double_t fraction = 0.;
      fMatch = 0;
      fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
      if(fSingleTrackEmb) {
	AliVParticle *vp = GetEmbeddedConstituent(jet1);
	if(vp) {
	  fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
	  fMatch = 1;
	}
      } else {
	jet2 = jet1->ClosestJet();
	fraction = jetCont->GetFractionSharedPt(jet1);
	fMatch = 1;
	if(fMinFractionShared>0.) {
	  if(fraction>fMinFractionShared) {
	    fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
	    fMatch = 1;
	  } else
	    fMatch = 0;
	}
      }

      //      if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();

      //Fill histograms for matched jets
      fh2MSubMatch[fCentBin]->Fill(var,fMatch);
      if(fMatch==1) {
	Double_t drToLJ = -1.;
	if(jetL) drToLJ = jet1->DeltaR(jetL);
	fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
	if(jet2) {
	  Double_t var2 = jet2->M();
	  if(fJetMassVarType==kRatMPt) {
	    if(jet2->Pt()>0. || jet2->Pt()<0.) var2 = jet2->M()/jet2->Pt();
	  }
          fh3MSubPtTrueLeadPt[fCentBin]->Fill(var,jet2->Pt(),jet1->MaxTrackPt());
          fh3MTruePtTrueLeadPt[fCentBin]->Fill(var2,jet2->Pt(),jet1->MaxTrackPt());
          fh3PtTrueDeltaMLeadPt[fCentBin]->Fill(jet2->Pt(),var-var2,jet1->MaxTrackPt());
          if(jet2->M()>0.) fh3PtTrueDeltaMRelLeadPt[fCentBin]->Fill(jet2->Pt(),(var-var2)/var2,jet1->MaxTrackPt());
	  Double_t varsp[5] = {var,var2,ptjet1,jet2->Pt(),jet1->MaxTrackPt()};
	  fhnMassResponse[fCentBin]->Fill(varsp);
	}
      }
      
      if(fCreateTree) {      
	fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
	if(jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
	else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
	fArea = (Float_t)jet1->Area();
	fAreaPhi = (Float_t)jet1->AreaPhi();
	fAreaEta = (Float_t)jet1->AreaEta();
	fRho  = (Float_t)jetCont->GetRhoVal();
	fRhoM = (Float_t)jetCont->GetRhoMassVal();
	fNConst = (Int_t)jet1->GetNumberOfTracks();
	fTreeJetBkg->Fill();
      }
    } //jet1 loop
  }//jetCont


  return kTRUE;
}

//________________________________________________________________________
AliVParticle* AliAnalysisTaskJetShapeConst::GetEmbeddedConstituent(AliEmcalJet *jet) {

  AliJetContainer *jetCont = GetJetContainer(fContainerBase);
  AliVParticle *vp = 0x0;
  AliVParticle *vpe = 0x0; //embedded particle
  Int_t nc = 0;
  for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
    vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
    if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
    if(!vpe) vpe = vp;
    else if(vp->Pt()>vpe->Pt()) vpe = vp;
    nc++;
  }
  AliDebug(11,Form("Found %d embedded particles",nc));
  return vpe;
}


//________________________________________________________________________
Bool_t AliAnalysisTaskJetShapeConst::RetrieveEventObjects() {
  //
  // retrieve event objects
  //

  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
    return kFALSE;

  AliJetContainer *jetCont = GetJetContainer(fContainerBase);
  jetCont->LoadRhoMass(InputEvent());

  return kTRUE;
}

//_______________________________________________________________________
void AliAnalysisTaskJetShapeConst::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
f4b02da3	1	//
	2	// Do subtraction for jet shapes using constituents arXiv:1403.3108
	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 <TF1.h>
	12	#include <TList.h>
	13	#include <TLorentzVector.h>
	14	#include <TProfile.h>
	15	#include <TChain.h>
	16	#include <TSystem.h>
	17	#include <TFile.h>
	18	#include <TKey.h>
	19	#include <TTree.h>
	20
	21	#include "AliVCluster.h"
	22	#include "AliVTrack.h"
	23	#include "AliEmcalJet.h"
	24	#include "AliRhoParameter.h"
	25	#include "AliLog.h"
	26	#include "AliEmcalParticle.h"
	27	#include "AliMCEvent.h"
	28	#include "AliAODEvent.h"
	29	#include "AliGenPythiaEventHeader.h"
	30	#include "AliAODMCHeader.h"
	31	#include "AliAnalysisManager.h"
	32	#include "AliJetContainer.h"
	33	#include "AliParticleContainer.h"
	34
	35	#include "AliAnalysisTaskJetShapeConst.h"
	36
	37	ClassImp(AliAnalysisTaskJetShapeConst)
	38
	39	//________________________________________________________________________
	40	AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst() :
	41	AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeConst", kTRUE),
	42	fContainerBase(0),
	43	fContainerSub(1),
a5a1c51a	44	fContainerNoEmb(2),
f4b02da3	45	fMinFractionShared(0),
	46	fSingleTrackEmb(kFALSE),
	47	fCreateTree(kFALSE),
d3660d18	48	fJetMassVarType(kMass),
80806b22	49	fUseSumw2(0),
f4b02da3	50	fTreeJetBkg(),
	51	fJet1Vec(new TLorentzVector()),
	52	fJet2Vec(new TLorentzVector()),
	53	fJetSubVec(new TLorentzVector()),
	54	fArea(0),
	55	fAreaPhi(0),
	56	fAreaEta(0),
	57	fRho(0),
	58	fRhoM(0),
	59	fNConst(0),
	60	fMatch(0),
	61	fh2MSubMatch(0x0),
	62	fh2MSubPtRawAll(0x0),
b09cf300	63	fh3MSubPtRawDRMatch(0x0),
dc481e30	64	fh3MSubPtTrueLeadPt(0x0),
	65	fh3MTruePtTrueLeadPt(0x0),
	66	fh3PtTrueDeltaMLeadPt(0x0),
	67	fh3PtTrueDeltaMRelLeadPt(0x0),
f4b02da3	68	fhnMassResponse(0x0)
	69	{
	70	// Default constructor.
	71
dc481e30	72	fh2MSubMatch = new TH2F*[fNcentBins];
	73	fh2MSubPtRawAll = new TH2F*[fNcentBins];
	74	fh3MSubPtRawDRMatch = new TH3F*[fNcentBins];
	75	fh3MSubPtTrueLeadPt = new TH3F*[fNcentBins];
	76	fh3MTruePtTrueLeadPt = new TH3F*[fNcentBins];
	77	fh3PtTrueDeltaMLeadPt = new TH3F*[fNcentBins];
	78	fh3PtTrueDeltaMRelLeadPt = new TH3F*[fNcentBins];
	79	fhnMassResponse = new THnSparse*[fNcentBins];
f4b02da3	80
f4b02da3	81	for (Int_t i = 0; i < fNcentBins; i++) {
dc481e30	82	fh2MSubMatch[i] = 0;
	83	fh2MSubPtRawAll[i] = 0;
	84	fh3MSubPtRawDRMatch[i] = 0;
	85	fh3MSubPtTrueLeadPt[i] = 0;
	86	fh3MTruePtTrueLeadPt[i] = 0;
	87	fh3PtTrueDeltaMLeadPt[i] = 0;
	88	fh3PtTrueDeltaMRelLeadPt[i] = 0;
	89	fhnMassResponse[i] = 0;
f4b02da3	90	}
	91
	92	SetMakeGeneralHistograms(kTRUE);
713c5683	93	if(fCreateTree) DefineOutput(2, TTree::Class());
f4b02da3	94	}
	95
	96	//________________________________________________________________________
	97	AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst(const char *name) :
	98	AliAnalysisTaskEmcalJet(name, kTRUE),
	99	fContainerBase(0),
	100	fContainerSub(1),
a5a1c51a	101	fContainerNoEmb(2),
f4b02da3	102	fMinFractionShared(0),
	103	fSingleTrackEmb(kFALSE),
	104	fCreateTree(kFALSE),
d3660d18	105	fJetMassVarType(kMass),
80806b22	106	fUseSumw2(0),
f4b02da3	107	fTreeJetBkg(0),
	108	fJet1Vec(new TLorentzVector()),
	109	fJet2Vec(new TLorentzVector()),
	110	fJetSubVec(new TLorentzVector()),
	111	fArea(0),
	112	fAreaPhi(0),
	113	fAreaEta(0),
	114	fRho(0),
	115	fRhoM(0),
	116	fNConst(0),
	117	fMatch(0),
	118	fh2MSubMatch(0x0),
	119	fh2MSubPtRawAll(0x0),
b09cf300	120	fh3MSubPtRawDRMatch(0x0),
dc481e30	121	fh3MSubPtTrueLeadPt(0x0),
	122	fh3MTruePtTrueLeadPt(0x0),
	123	fh3PtTrueDeltaMLeadPt(0x0),
	124	fh3PtTrueDeltaMRelLeadPt(0x0),
f4b02da3	125	fhnMassResponse(0x0)
	126	{
	127	// Standard constructor.
	128
dc481e30	129	fh2MSubMatch = new TH2F*[fNcentBins];
	130	fh2MSubPtRawAll = new TH2F*[fNcentBins];
	131	fh3MSubPtRawDRMatch = new TH3F*[fNcentBins];
	132	fh3MSubPtTrueLeadPt = new TH3F*[fNcentBins];
	133	fh3MTruePtTrueLeadPt = new TH3F*[fNcentBins];
	134	fh3PtTrueDeltaMLeadPt = new TH3F*[fNcentBins];
	135	fh3PtTrueDeltaMRelLeadPt = new TH3F*[fNcentBins];
	136	fhnMassResponse = new THnSparse*[fNcentBins];
f4b02da3	137
f4b02da3	138	for (Int_t i = 0; i < fNcentBins; i++) {
dc481e30	139	fh2MSubMatch[i] = 0;
	140	fh2MSubPtRawAll[i] = 0;
	141	fh3MSubPtRawDRMatch[i] = 0;
	142	fh3MSubPtTrueLeadPt[i] = 0;
	143	fh3MTruePtTrueLeadPt[i] = 0;
	144	fh3PtTrueDeltaMLeadPt[i] = 0;
	145	fh3PtTrueDeltaMRelLeadPt[i] = 0;
	146	fhnMassResponse[i] = 0;
f4b02da3	147	}
	148
	149	SetMakeGeneralHistograms(kTRUE);
713c5683	150	if(fCreateTree) DefineOutput(2, TTree::Class());
f4b02da3	151	}
	152
	153	//________________________________________________________________________
	154	AliAnalysisTaskJetShapeConst::~AliAnalysisTaskJetShapeConst()
	155	{
	156	// Destructor.
	157	}
	158
	159	//________________________________________________________________________
	160	void AliAnalysisTaskJetShapeConst::UserCreateOutputObjects()
	161	{
	162	// Create user output.
	163
	164	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
	165
	166	Bool_t oldStatus = TH1::AddDirectoryStatus();
	167	TH1::AddDirectory(kFALSE);
	168
	169	const Int_t nBinsPt = 200;
	170	const Double_t minPt = -50.;
	171	const Double_t maxPt = 150.;
	172
d3660d18	173	Int_t nBinsM = 100;
	174	Double_t minM = -20.;
	175	Double_t maxM = 80.;
	176	if(fJetMassVarType==kRatMPt) {
	177	nBinsM = 100;
	178	minM = -0.2;
	179	maxM = 0.8;
	180	}
f4b02da3	181
d3660d18	182	Int_t nBinsDM = 100;
a3c67475	183	Double_t minDM = -25.;
a3c67475	184	Double_t maxDM = 25.;
d3660d18	185	if(fJetMassVarType==kRatMPt) {
	186	nBinsDM = 100;
	187	minDM = -0.5;
	188	maxDM = 0.5;
	189	}
a5a1c51a	190
	191	const Int_t nBinsDRToLJ = 20; //distance to leading jet in Pb-Pb only event
	192	const Double_t minDRToLJ = 0.;
	193	const Double_t maxDRToLJ = 1.;
8d95aff1	194
dc481e30	195	const Int_t nBinsPtLead = 20;
	196	const Double_t minPtLead = 0.;
	197	const Double_t maxPtLead = 20.;
	198
f4b02da3	199	//Binning for THnSparse
8d95aff1	200	const Int_t nBinsSparse0 = 5;
dc481e30	201	const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsPtLead};
	202	const Double_t xmin0[nBinsSparse0] = { minM, minM, minPt, minPt, minPtLead};
	203	const Double_t xmax0[nBinsSparse0] = { maxM, maxM, maxPt, maxPt, maxPtLead};
f4b02da3	204
	205	TString histName = "";
	206	TString histTitle = "";
d3660d18	207	TString varName = "#it{M}_{jet}";
	208	if(fJetMassVarType==kRatMPt) varName = "#it{M}_{jet}/#it{p}_{T,jet}";
	209
f4b02da3	210	for (Int_t i = 0; i < fNcentBins; i++) {
f4b02da3	211	histName = Form("fh2MSubMatch_%d",i);
d3660d18	212	histTitle = Form("fh2MSubMatch_%d;%s;match",i,varName.Data());
f4b02da3	213	fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
	214	fOutput->Add(fh2MSubMatch[i]);
	215
	216	histName = Form("fh2MSubPtRawAll_%d",i);
d3660d18	217	histTitle = Form("fh2MSubPtRawAll_%d;%s;#it{p}_{T}",i,varName.Data());
f4b02da3	218	fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
	219	fOutput->Add(fh2MSubPtRawAll[i]);
	220
b09cf300	221	histName = Form("fh3MSubPtRawDRMatch_%d",i);
d3660d18	222	histTitle = Form("fh3MSubPtRawDRMatch_%d;%s;#it{p}_{T}",i,varName.Data());
b09cf300	223	fh3MSubPtRawDRMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
b09cf300	224	fOutput->Add(fh3MSubPtRawDRMatch[i]);
f4b02da3	225
dc481e30	226	histName = Form("fh3MSubPtTrueLeadPt_%d",i);
	227	histTitle = Form("fh3MSubPtTrueLeadPt_%d;%s;#it{p}_{T}",i,varName.Data());
	228	fh3MSubPtTrueLeadPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsPtLead,minPtLead,maxPtLead);
	229	fOutput->Add(fh3MSubPtTrueLeadPt[i]);
f4b02da3	230
dc481e30	231	histName = Form("fh3MTruePtTrueLeadPt_%d",i);
	232	histTitle = Form("fh3MTruePtTrueLeadPt_%d;%s;#it{p}_{T}",i,varName.Data());
	233	fh3MTruePtTrueLeadPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsPtLead,minPtLead,maxPtLead);
	234	fOutput->Add(fh3MTruePtTrueLeadPt[i]);
f4b02da3	235
dc481e30	236	histName = Form("fh3PtTrueDeltaMLeadPt_%d",i);
	237	histTitle = Form("fh3PtTrueDeltaMLeadPt_%d;#it{p}_{T,true};#Delta %s",i,varName.Data());
	238	fh3PtTrueDeltaMLeadPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDM,minDM,maxDM,nBinsPtLead,minPtLead,maxPtLead);
	239	fOutput->Add(fh3PtTrueDeltaMLeadPt[i]);
f4b02da3	240
dc481e30	241	histName = Form("fh3PtTrueDeltaMRelLeadPt_%d",i);
	242	histTitle = Form("fh3PtTrueDeltaMRelLeadPt_%d;#it{p}_{T,true};Rel #Delta %s",i,varName.Data());
	243	fh3PtTrueDeltaMRelLeadPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,400,-1.,3.,nBinsPtLead,minPtLead,maxPtLead);
	244	fOutput->Add(fh3PtTrueDeltaMRelLeadPt[i]);
0947b103	245
f4b02da3	246	histName = Form("fhnMassResponse_%d",i);
dc481e30	247	histTitle = Form("fhnMassResponse_%d;%s sub;%s true;#it{p}_{T,sub};#it{p}_{T,true};#it{p}_{T,lead trk}",i,varName.Data(),varName.Data());
f4b02da3	248	fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
	249	fOutput->Add(fhnMassResponse[i]);
	250	}
	251
80806b22	252	if(fUseSumw2) {
	253	// =========== Switch on Sumw2 for all histos ===========
	254	for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
	255	TH1 h1 = dynamic_cast<TH1>(fOutput->At(i));
	256	if (h1){
	257	h1->Sumw2();
	258	continue;
	259	}
	260	THnSparse hn = dynamic_cast<THnSparse>(fOutput->At(i));
	261	if(hn)hn->Sumw2();
f4b02da3	262	}
f4b02da3	263	}
	264
	265	TH1::AddDirectory(oldStatus);
	266
	267	// Create a tree.
	268	if(fCreateTree) {
	269	fTreeJetBkg = new TTree("fTreeJetSubConst", "fTreeJetSubConst");
	270	fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
	271	fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
	272	fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
	273	fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
	274	fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
	275	fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
	276	fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
	277	fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
	278	fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
	279	}
	280
	281	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
	282	if(fCreateTree) PostData(2, fTreeJetBkg);
	283	}
	284
	285	//________________________________________________________________________
	286	Bool_t AliAnalysisTaskJetShapeConst::Run()
	287	{
	288	// Run analysis code here, if needed. It will be executed before FillHistograms().
	289
	290	return kTRUE;
	291	}
	292
	293	//________________________________________________________________________
	294	Bool_t AliAnalysisTaskJetShapeConst::FillHistograms()
	295	{
	296	// Fill histograms.
	297
8d95aff1	298	AliEmcalJet* jet1 = NULL; //AA jet
8d95aff1	299	AliEmcalJet *jet2 = NULL; //Embedded Pythia jet
a5a1c51a	300	// AliEmcalJet *jet1T = NULL; //tagged AA jet
8d95aff1	301	// AliEmcalJet *jet2T = NULL; //tagged Pythia jet
8d95aff1	302	AliEmcalJet *jetS = NULL; //subtracted jet
f4b02da3	303	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
f4b02da3	304	AliJetContainer *jetContS = GetJetContainer(fContainerSub);
a5a1c51a	305
	306	//Get leading jet in Pb-Pb event without embedded objects
	307	AliJetContainer *jetContNoEmb = GetJetContainer(fContainerNoEmb);
	308	AliEmcalJet *jetL = NULL;
	309	if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");
	310
f4b02da3	311	AliDebug(11,Form("NJets Incl: %d Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
	312	if(jetCont) {
	313	jetCont->ResetCurrentID();
	314	while((jet1 = jetCont->GetNextAcceptJet())) {
8d95aff1	315	jet2 = NULL;
8d95aff1	316	jetS = NULL;
d3660d18	317
f4b02da3	318	//Get constituent subtacted version of jet
	319	Int_t ifound = 0;
	320	Int_t ilab = -1;
	321	for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
	322	//if(ifound==1) continue;
	323	jetS = jetContS->GetJet(i);
	324	if(jetS->GetLabel()==jet1->GetLabel()) { // && jetS->Pt()>0.) {
	325	ifound++;
	326	if(ifound==1) ilab = i;
	327	}
	328	}
	329	if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
	330	if(ifound==0) jetS = 0x0;
	331	else jetS = jetContS->GetJet(ilab);
	332
d3660d18	333	Double_t mjetS = jetS->M();
	334	Double_t ptjet1 = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
	335	Double_t var = mjetS;
	336	if(fJetMassVarType==kRatMPt) {
	337	if(ptjet1>0. \|\| ptjet1<0.) var = mjetS/ptjet1;
	338	else var = -999.;
	339	}
	340
f4b02da3	341	//Fill histograms for all AA jets
d3660d18	342	fh2MSubPtRawAll[fCentBin]->Fill(var,ptjet1);
f4b02da3	343
	344	Double_t fraction = 0.;
	345	fMatch = 0;
	346	fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
	347	if(fSingleTrackEmb) {
	348	AliVParticle *vp = GetEmbeddedConstituent(jet1);
	349	if(vp) {
	350	fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
	351	fMatch = 1;
	352	}
	353	} else {
	354	jet2 = jet1->ClosestJet();
	355	fraction = jetCont->GetFractionSharedPt(jet1);
	356	fMatch = 1;
	357	if(fMinFractionShared>0.) {
	358	if(fraction>fMinFractionShared) {
	359	fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
	360	fMatch = 1;
	361	} else
	362	fMatch = 0;
	363	}
	364	}
	365
8d95aff1	366	// if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
8d95aff1	367
f4b02da3	368	//Fill histograms for matched jets
d3660d18	369	fh2MSubMatch[fCentBin]->Fill(var,fMatch);
f4b02da3	370	if(fMatch==1) {
b09cf300	371	Double_t drToLJ = -1.;
b09cf300	372	if(jetL) drToLJ = jet1->DeltaR(jetL);
d3660d18	373	fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
f4b02da3	374	if(jet2) {
e0036e47	375	Double_t var2 = jet2->M();
	376	if(fJetMassVarType==kRatMPt) {
	377	if(jet2->Pt()>0. \|\| jet2->Pt()<0.) var2 = jet2->M()/jet2->Pt();
	378	}
dc481e30	379	fh3MSubPtTrueLeadPt[fCentBin]->Fill(var,jet2->Pt(),jet1->MaxTrackPt());
	380	fh3MTruePtTrueLeadPt[fCentBin]->Fill(var2,jet2->Pt(),jet1->MaxTrackPt());
	381	fh3PtTrueDeltaMLeadPt[fCentBin]->Fill(jet2->Pt(),var-var2,jet1->MaxTrackPt());
	382	if(jet2->M()>0.) fh3PtTrueDeltaMRelLeadPt[fCentBin]->Fill(jet2->Pt(),(var-var2)/var2,jet1->MaxTrackPt());
	383	Double_t varsp[5] = {var,var2,ptjet1,jet2->Pt(),jet1->MaxTrackPt()};
d3660d18	384	fhnMassResponse[fCentBin]->Fill(varsp);
f4b02da3	385	}
	386	}
	387
	388	if(fCreateTree) {
	389	fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
c090d64e	390	if(jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
f4b02da3	391	else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
	392	fArea = (Float_t)jet1->Area();
	393	fAreaPhi = (Float_t)jet1->AreaPhi();
	394	fAreaEta = (Float_t)jet1->AreaEta();
	395	fRho = (Float_t)jetCont->GetRhoVal();
	396	fRhoM = (Float_t)jetCont->GetRhoMassVal();
	397	fNConst = (Int_t)jet1->GetNumberOfTracks();
	398	fTreeJetBkg->Fill();
	399	}
	400	} //jet1 loop
	401	}//jetCont
	402
	403
	404	return kTRUE;
	405	}
	406
	407	//________________________________________________________________________
	408	AliVParticle* AliAnalysisTaskJetShapeConst::GetEmbeddedConstituent(AliEmcalJet *jet) {
	409
	410	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
	411	AliVParticle *vp = 0x0;
	412	AliVParticle *vpe = 0x0; //embedded particle
	413	Int_t nc = 0;
	414	for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
	415	vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
	416	if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
	417	if(!vpe) vpe = vp;
	418	else if(vp->Pt()>vpe->Pt()) vpe = vp;
	419	nc++;
	420	}
	421	AliDebug(11,Form("Found %d embedded particles",nc));
	422	return vpe;
	423	}
	424
	425
	426	//________________________________________________________________________
	427	Bool_t AliAnalysisTaskJetShapeConst::RetrieveEventObjects() {
	428	//
	429	// retrieve event objects
	430	//
	431
	432	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
	433	return kFALSE;
	434
	435	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
	436	jetCont->LoadRhoMass(InputEvent());
	437
	438	return kTRUE;
	439	}
	440
	441	//_______________________________________________________________________
	442	void AliAnalysisTaskJetShapeConst::Terminate(Option_t *)
	443	{
	444	// Called once at the end of the analysis.
	445	}
	446