[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),
  fMinFractionShared(0),
  fSingleTrackEmb(kFALSE),
  fCreateTree(kFALSE),
  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),
  fh2MSubPtRawMatch(0x0),
  fh2MSubPtTrue(0x0),
  fh2MTruePtTrue(0x0),
  fh2PtTrueDeltaM(0x0),
  fh2PtTrueDeltaMRel(0x0),
  fhnMassResponse(0x0)
{
  // Default constructor.

  fh2MSubMatch        = new TH2F*[fNcentBins];
  fh2MSubPtRawAll     = new TH2F*[fNcentBins];
  fh2MSubPtRawMatch   = new TH2F*[fNcentBins];
  fh2MSubPtTrue       = new TH2F*[fNcentBins];
  fh2MTruePtTrue      = new TH2F*[fNcentBins];
  fh2PtTrueDeltaM     = new TH2F*[fNcentBins];
  fh2PtTrueDeltaMRel  = new TH2F*[fNcentBins];
  fhnMassResponse     = new THnSparse*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh2MSubMatch[i]        = 0;
    fh2MSubPtRawAll[i]     = 0;
    fh2MSubPtRawMatch[i]   = 0;
    fh2MSubPtTrue[i]       = 0;
    fh2MTruePtTrue[i]      = 0;
    fh2PtTrueDeltaM[i]     = 0;
    fh2PtTrueDeltaMRel[i]  = 0;
    fhnMassResponse[i]     = 0;
  }

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

//________________________________________________________________________
AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),  
  fContainerBase(0),
  fContainerSub(1),
  fMinFractionShared(0),
  fSingleTrackEmb(kFALSE),
  fCreateTree(kFALSE),
  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),
  fh2MSubPtRawMatch(0x0),
  fh2MSubPtTrue(0x0),
  fh2MTruePtTrue(0x0),
  fh2PtTrueDeltaM(0x0),
  fh2PtTrueDeltaMRel(0x0),
  fhnMassResponse(0x0)
{
  // Standard constructor.

  fh2MSubMatch        = new TH2F*[fNcentBins];
  fh2MSubPtRawAll     = new TH2F*[fNcentBins];
  fh2MSubPtRawMatch   = new TH2F*[fNcentBins];
  fh2MSubPtTrue       = new TH2F*[fNcentBins];
  fh2MTruePtTrue      = new TH2F*[fNcentBins];
  fh2PtTrueDeltaM     = new TH2F*[fNcentBins];
  fh2PtTrueDeltaMRel  = new TH2F*[fNcentBins];
  fhnMassResponse     = new THnSparse*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh2MSubMatch[i]        = 0;
    fh2MSubPtRawAll[i]     = 0;
    fh2MSubPtRawMatch[i]   = 0;
    fh2MSubPtTrue[i]       = 0;
    fh2MTruePtTrue[i]      = 0;
    fh2PtTrueDeltaM[i]     = 0;
    fh2PtTrueDeltaMRel[i]  = 0;
    fhnMassResponse[i]     = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
  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.;

  const Int_t nBinsM  = 150;
  const Double_t minM = -50.;
  const Double_t maxM = 100.;

  const Int_t nBinsMT  = 50;
  const Double_t minMT = 0.;
  const Double_t maxMT = 50.;

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

  TString histName = "";
  TString histTitle = "";
  for (Int_t i = 0; i < fNcentBins; i++) {
    histName = Form("fh2MSubMatch_%d",i);
    histTitle = Form("fh2MSubMatch_%d;#it{M}_{sub};match",i);
    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;#it{M}_{sub};#it{p}_{T}",i);
    fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
    fOutput->Add(fh2MSubPtRawAll[i]);

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

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

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

    histName = Form("fh2PtTrueDeltaM_%d",i);
    histTitle = Form("fh2PtTrueDeltaM_%d;#it{p}_{T,true};#it{M}_{sub}-#it{M}_{true}",i);
    fh2PtTrueDeltaM[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,-50.,50.);
    fOutput->Add(fh2PtTrueDeltaM[i]);

    histName = Form("fh2PtTrueDeltaMRel_%d",i);
    histTitle = Form("fh2PtTrueDeltaMRel_%d;#it{p}_{T,true};(#it{M}_{sub}-#it{M}_{true})/#it{M}_{true}",i);
    fh2PtTrueDeltaMRel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,200,-1.,1.);
    fOutput->Add(fh2PtTrueDeltaMRel[i]);

    histName = Form("fhnMassResponse_%d",i);
    histTitle = Form("fhnMassResponse_%d;#it{M}_{sub};#it{M}_{true};#it{p}_{T,sub};#it{p}_{T,true};#it{M}_{sub}^{tagged}",i);
    fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
    fOutput->Add(fhnMassResponse[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);

  // 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);
  AliDebug(11,Form("NJets  Incl: %d  Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
  if(jetCont) {
    jetCont->ResetCurrentID();
    while((jet1 = jetCont->GetNextAcceptJet())) {
      jet2  = NULL;
      jet1T = NULL;
      //   jet2T = NULL;
      jetS  = NULL;
      if(jet1->GetTagStatus()>0) jet1T = jet1->GetTaggedJet();
 
      //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);

      //Fill histograms for all AA jets
      fh2MSubPtRawAll[fCentBin]->Fill(jetS->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area());

      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(jetS->M(),fMatch);
      if(fMatch==1) {
	fh2MSubPtRawMatch[fCentBin]->Fill(jetS->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area());
	if(jet2) {
	  fh2MSubPtTrue[fCentBin]->Fill(jetS->M(),jet2->Pt());
	  fh2MTruePtTrue[fCentBin]->Fill(jet2->M(),jet2->Pt());
	  fh2PtTrueDeltaM[fCentBin]->Fill(jet2->Pt(),jetS->M()-jet2->M());
	  if(jet2->M()>0.) fh2PtTrueDeltaMRel[fCentBin]->Fill(jet2->Pt(),(jetS->M()-jet2->M())/jet2->M());
	  Double_t mJet1Tagged = -1.;
	  if(jet1T) mJet1Tagged = jet1T->M();
	  Double_t var[5] = {jetS->M(),jet2->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area(),jet2->Pt(),mJet1Tagged};
	  fhnMassResponse[fCentBin]->Fill(var);
	}
      }
      
      if(fCreateTree) {      
	fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
	if(jetS && 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),
	44	fMinFractionShared(0),
	45	fSingleTrackEmb(kFALSE),
	46	fCreateTree(kFALSE),
	47	fTreeJetBkg(),
	48	fJet1Vec(new TLorentzVector()),
	49	fJet2Vec(new TLorentzVector()),
	50	fJetSubVec(new TLorentzVector()),
	51	fArea(0),
	52	fAreaPhi(0),
	53	fAreaEta(0),
	54	fRho(0),
	55	fRhoM(0),
	56	fNConst(0),
	57	fMatch(0),
	58	fh2MSubMatch(0x0),
	59	fh2MSubPtRawAll(0x0),
	60	fh2MSubPtRawMatch(0x0),
	61	fh2MSubPtTrue(0x0),
	62	fh2MTruePtTrue(0x0),
	63	fh2PtTrueDeltaM(0x0),
0947b103	64	fh2PtTrueDeltaMRel(0x0),
f4b02da3	65	fhnMassResponse(0x0)
	66	{
	67	// Default constructor.
	68
	69	fh2MSubMatch = new TH2F*[fNcentBins];
	70	fh2MSubPtRawAll = new TH2F*[fNcentBins];
	71	fh2MSubPtRawMatch = new TH2F*[fNcentBins];
	72	fh2MSubPtTrue = new TH2F*[fNcentBins];
	73	fh2MTruePtTrue = new TH2F*[fNcentBins];
	74	fh2PtTrueDeltaM = new TH2F*[fNcentBins];
0947b103	75	fh2PtTrueDeltaMRel = new TH2F*[fNcentBins];
f4b02da3	76	fhnMassResponse = new THnSparse*[fNcentBins];
	77
	78	for (Int_t i = 0; i < fNcentBins; i++) {
	79	fh2MSubMatch[i] = 0;
	80	fh2MSubPtRawAll[i] = 0;
	81	fh2MSubPtRawMatch[i] = 0;
	82	fh2MSubPtTrue[i] = 0;
	83	fh2MTruePtTrue[i] = 0;
	84	fh2PtTrueDeltaM[i] = 0;
0947b103	85	fh2PtTrueDeltaMRel[i] = 0;
f4b02da3	86	fhnMassResponse[i] = 0;
	87	}
	88
	89	SetMakeGeneralHistograms(kTRUE);
	90	DefineOutput(2, TTree::Class());
	91	}
	92
	93	//________________________________________________________________________
	94	AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst(const char *name) :
	95	AliAnalysisTaskEmcalJet(name, kTRUE),
	96	fContainerBase(0),
	97	fContainerSub(1),
	98	fMinFractionShared(0),
	99	fSingleTrackEmb(kFALSE),
	100	fCreateTree(kFALSE),
	101	fTreeJetBkg(0),
	102	fJet1Vec(new TLorentzVector()),
	103	fJet2Vec(new TLorentzVector()),
	104	fJetSubVec(new TLorentzVector()),
	105	fArea(0),
	106	fAreaPhi(0),
	107	fAreaEta(0),
	108	fRho(0),
	109	fRhoM(0),
	110	fNConst(0),
	111	fMatch(0),
	112	fh2MSubMatch(0x0),
	113	fh2MSubPtRawAll(0x0),
	114	fh2MSubPtRawMatch(0x0),
	115	fh2MSubPtTrue(0x0),
	116	fh2MTruePtTrue(0x0),
	117	fh2PtTrueDeltaM(0x0),
0947b103	118	fh2PtTrueDeltaMRel(0x0),
f4b02da3	119	fhnMassResponse(0x0)
	120	{
	121	// Standard constructor.
	122
	123	fh2MSubMatch = new TH2F*[fNcentBins];
	124	fh2MSubPtRawAll = new TH2F*[fNcentBins];
	125	fh2MSubPtRawMatch = new TH2F*[fNcentBins];
	126	fh2MSubPtTrue = new TH2F*[fNcentBins];
	127	fh2MTruePtTrue = new TH2F*[fNcentBins];
	128	fh2PtTrueDeltaM = new TH2F*[fNcentBins];
0947b103	129	fh2PtTrueDeltaMRel = new TH2F*[fNcentBins];
f4b02da3	130	fhnMassResponse = new THnSparse*[fNcentBins];
	131
	132	for (Int_t i = 0; i < fNcentBins; i++) {
	133	fh2MSubMatch[i] = 0;
	134	fh2MSubPtRawAll[i] = 0;
	135	fh2MSubPtRawMatch[i] = 0;
	136	fh2MSubPtTrue[i] = 0;
	137	fh2MTruePtTrue[i] = 0;
	138	fh2PtTrueDeltaM[i] = 0;
0947b103	139	fh2PtTrueDeltaMRel[i] = 0;
f4b02da3	140	fhnMassResponse[i] = 0;
	141	}
	142
	143	SetMakeGeneralHistograms(kTRUE);
	144	DefineOutput(2, TTree::Class());
	145	}
	146
	147	//________________________________________________________________________
	148	AliAnalysisTaskJetShapeConst::~AliAnalysisTaskJetShapeConst()
	149	{
	150	// Destructor.
	151	}
	152
	153	//________________________________________________________________________
	154	void AliAnalysisTaskJetShapeConst::UserCreateOutputObjects()
	155	{
	156	// Create user output.
	157
	158	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
	159
	160	Bool_t oldStatus = TH1::AddDirectoryStatus();
	161	TH1::AddDirectory(kFALSE);
	162
	163	const Int_t nBinsPt = 200;
	164	const Double_t minPt = -50.;
	165	const Double_t maxPt = 150.;
	166
	167	const Int_t nBinsM = 150;
	168	const Double_t minM = -50.;
	169	const Double_t maxM = 100.;
	170
8d95aff1	171	const Int_t nBinsMT = 50;
	172	const Double_t minMT = 0.;
	173	const Double_t maxMT = 50.;
	174
f4b02da3	175	//Binning for THnSparse
8d95aff1	176	const Int_t nBinsSparse0 = 5;
	177	const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsMT};
	178	const Double_t xmin0[nBinsSparse0] = { minM, minM, minPt, minPt, minMT};
	179	const Double_t xmax0[nBinsSparse0] = { maxM, maxM, maxPt, maxPt, maxMT};
f4b02da3	180
	181	TString histName = "";
	182	TString histTitle = "";
	183	for (Int_t i = 0; i < fNcentBins; i++) {
	184	histName = Form("fh2MSubMatch_%d",i);
	185	histTitle = Form("fh2MSubMatch_%d;#it{M}_{sub};match",i);
	186	fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
	187	fOutput->Add(fh2MSubMatch[i]);
	188
	189	histName = Form("fh2MSubPtRawAll_%d",i);
	190	histTitle = Form("fh2MSubPtRawAll_%d;#it{M}_{sub};#it{p}_{T}",i);
	191	fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
	192	fOutput->Add(fh2MSubPtRawAll[i]);
	193
	194	histName = Form("fh2MSubPtRawMatch_%d",i);
	195	histTitle = Form("fh2MSubPtRawMatch_%d;#it{M}_{sub};#it{p}_{T}",i);
	196	fh2MSubPtRawMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
	197	fOutput->Add(fh2MSubPtRawMatch[i]);
	198
	199	histName = Form("fh2MSubPtTrue_%d",i);
	200	histTitle = Form("fh2MSubPtTrue_%d;#it{M}_{sub};#it{p}_{T}",i);
	201	fh2MSubPtTrue[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
	202	fOutput->Add(fh2MSubPtTrue[i]);
	203
	204	histName = Form("fh2MTruePtTrue_%d",i);
	205	histTitle = Form("fh2MTruePtTrue_%d;#it{M}_{sub};#it{p}_{T}",i);
	206	fh2MTruePtTrue[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
	207	fOutput->Add(fh2MTruePtTrue[i]);
	208
	209	histName = Form("fh2PtTrueDeltaM_%d",i);
	210	histTitle = Form("fh2PtTrueDeltaM_%d;#it{p}_{T,true};#it{M}_{sub}-#it{M}_{true}",i);
	211	fh2PtTrueDeltaM[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,-50.,50.);
	212	fOutput->Add(fh2PtTrueDeltaM[i]);
	213
0947b103	214	histName = Form("fh2PtTrueDeltaMRel_%d",i);
	215	histTitle = Form("fh2PtTrueDeltaMRel_%d;#it{p}_{T,true};(#it{M}_{sub}-#it{M}_{true})/#it{M}_{true}",i);
	216	fh2PtTrueDeltaMRel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,200,-1.,1.);
	217	fOutput->Add(fh2PtTrueDeltaMRel[i]);
	218
f4b02da3	219	histName = Form("fhnMassResponse_%d",i);
8d95aff1	220	histTitle = Form("fhnMassResponse_%d;#it{M}_{sub};#it{M}_{true};#it{p}_{T,sub};#it{p}_{T,true};#it{M}_{sub}^{tagged}",i);
f4b02da3	221	fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
	222	fOutput->Add(fhnMassResponse[i]);
	223	}
	224
	225	// =========== Switch on Sumw2 for all histos ===========
	226	for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
	227	TH1 h1 = dynamic_cast<TH1>(fOutput->At(i));
	228	if (h1){
	229	h1->Sumw2();
	230	continue;
	231	}
	232	THnSparse hn = dynamic_cast<THnSparse>(fOutput->At(i));
	233	if(hn)hn->Sumw2();
	234	}
	235
	236	TH1::AddDirectory(oldStatus);
	237
	238	// Create a tree.
	239	if(fCreateTree) {
	240	fTreeJetBkg = new TTree("fTreeJetSubConst", "fTreeJetSubConst");
	241	fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
	242	fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
	243	fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
	244	fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
	245	fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
	246	fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
	247	fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
	248	fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
	249	fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
	250	}
	251
	252	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
	253	if(fCreateTree) PostData(2, fTreeJetBkg);
	254	}
	255
	256	//________________________________________________________________________
	257	Bool_t AliAnalysisTaskJetShapeConst::Run()
	258	{
	259	// Run analysis code here, if needed. It will be executed before FillHistograms().
	260
	261	return kTRUE;
	262	}
	263
	264	//________________________________________________________________________
	265	Bool_t AliAnalysisTaskJetShapeConst::FillHistograms()
	266	{
	267	// Fill histograms.
	268
8d95aff1	269	AliEmcalJet* jet1 = NULL; //AA jet
	270	AliEmcalJet *jet2 = NULL; //Embedded Pythia jet
	271	AliEmcalJet *jet1T = NULL; //tagged AA jet
	272	// AliEmcalJet *jet2T = NULL; //tagged Pythia jet
	273	AliEmcalJet *jetS = NULL; //subtracted jet
f4b02da3	274	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
	275	AliJetContainer *jetContS = GetJetContainer(fContainerSub);
	276	AliDebug(11,Form("NJets Incl: %d Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
	277	if(jetCont) {
	278	jetCont->ResetCurrentID();
	279	while((jet1 = jetCont->GetNextAcceptJet())) {
8d95aff1	280	jet2 = NULL;
	281	jet1T = NULL;
	282	// jet2T = NULL;
	283	jetS = NULL;
	284	if(jet1->GetTagStatus()>0) jet1T = jet1->GetTaggedJet();
	285
f4b02da3	286	//Get constituent subtacted version of jet
	287	Int_t ifound = 0;
	288	Int_t ilab = -1;
	289	for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
	290	//if(ifound==1) continue;
	291	jetS = jetContS->GetJet(i);
	292	if(jetS->GetLabel()==jet1->GetLabel()) { // && jetS->Pt()>0.) {
	293	ifound++;
	294	if(ifound==1) ilab = i;
	295	}
	296	}
	297	if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
	298	if(ifound==0) jetS = 0x0;
	299	else jetS = jetContS->GetJet(ilab);
	300
	301	//Fill histograms for all AA jets
	302	fh2MSubPtRawAll[fCentBin]->Fill(jetS->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area());
	303
	304	Double_t fraction = 0.;
	305	fMatch = 0;
	306	fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
	307	if(fSingleTrackEmb) {
	308	AliVParticle *vp = GetEmbeddedConstituent(jet1);
	309	if(vp) {
	310	fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
	311	fMatch = 1;
	312	}
	313	} else {
	314	jet2 = jet1->ClosestJet();
	315	fraction = jetCont->GetFractionSharedPt(jet1);
	316	fMatch = 1;
	317	if(fMinFractionShared>0.) {
	318	if(fraction>fMinFractionShared) {
	319	fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
	320	fMatch = 1;
	321	} else
	322	fMatch = 0;
	323	}
	324	}
	325
8d95aff1	326	// if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
8d95aff1	327
f4b02da3	328	//Fill histograms for matched jets
	329	fh2MSubMatch[fCentBin]->Fill(jetS->M(),fMatch);
	330	if(fMatch==1) {
	331	fh2MSubPtRawMatch[fCentBin]->Fill(jetS->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area());
	332	if(jet2) {
	333	fh2MSubPtTrue[fCentBin]->Fill(jetS->M(),jet2->Pt());
	334	fh2MTruePtTrue[fCentBin]->Fill(jet2->M(),jet2->Pt());
	335	fh2PtTrueDeltaM[fCentBin]->Fill(jet2->Pt(),jetS->M()-jet2->M());
0947b103	336	if(jet2->M()>0.) fh2PtTrueDeltaMRel[fCentBin]->Fill(jet2->Pt(),(jetS->M()-jet2->M())/jet2->M());
8d95aff1	337	Double_t mJet1Tagged = -1.;
	338	if(jet1T) mJet1Tagged = jet1T->M();
	339	Double_t var[5] = {jetS->M(),jet2->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area(),jet2->Pt(),mJet1Tagged};
f4b02da3	340	fhnMassResponse[fCentBin]->Fill(var);
	341	}
	342	}
	343
	344	if(fCreateTree) {
	345	fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
	346	if(jetS && jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
	347	else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
	348	fArea = (Float_t)jet1->Area();
	349	fAreaPhi = (Float_t)jet1->AreaPhi();
	350	fAreaEta = (Float_t)jet1->AreaEta();
	351	fRho = (Float_t)jetCont->GetRhoVal();
	352	fRhoM = (Float_t)jetCont->GetRhoMassVal();
	353	fNConst = (Int_t)jet1->GetNumberOfTracks();
	354	fTreeJetBkg->Fill();
	355	}
	356	} //jet1 loop
	357	}//jetCont
	358
	359
	360	return kTRUE;
	361	}
	362
	363	//________________________________________________________________________
	364	AliVParticle* AliAnalysisTaskJetShapeConst::GetEmbeddedConstituent(AliEmcalJet *jet) {
	365
	366	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
	367	AliVParticle *vp = 0x0;
	368	AliVParticle *vpe = 0x0; //embedded particle
	369	Int_t nc = 0;
	370	for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
	371	vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
	372	if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
	373	if(!vpe) vpe = vp;
	374	else if(vp->Pt()>vpe->Pt()) vpe = vp;
	375	nc++;
	376	}
	377	AliDebug(11,Form("Found %d embedded particles",nc));
	378	return vpe;
	379	}
	380
	381
	382	//________________________________________________________________________
	383	Bool_t AliAnalysisTaskJetShapeConst::RetrieveEventObjects() {
	384	//
	385	// retrieve event objects
	386	//
	387
	388	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
	389	return kFALSE;
	390
	391	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
	392	jetCont->LoadRhoMass(InputEvent());
	393
	394	return kTRUE;
	395	}
	396
	397	//_______________________________________________________________________
	398	void AliAnalysisTaskJetShapeConst::Terminate(Option_t *)
	399	{
	400	// Called once at the end of the analysis.
	401	}
	402