[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),
  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),
  fh3MSubPtTrueDR(0x0),
  fh3MTruePtTrueDR(0x0),
  fh3PtTrueDeltaMDR(0x0),
  fh3PtTrueDeltaMRelDR(0x0),
  fhnMassResponse(0x0)
{
  // Default constructor.

  fh2MSubMatch         = new TH2F*[fNcentBins];
  fh2MSubPtRawAll      = new TH2F*[fNcentBins];
  fh3MSubPtRawDRMatch  = new TH3F*[fNcentBins];
  fh3MSubPtTrueDR      = new TH3F*[fNcentBins];
  fh3MTruePtTrueDR     = new TH3F*[fNcentBins];
  fh3PtTrueDeltaMDR    = new TH3F*[fNcentBins];
  fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
  fhnMassResponse      = new THnSparse*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh2MSubMatch[i]         = 0;
    fh2MSubPtRawAll[i]      = 0;
    fh3MSubPtRawDRMatch[i]  = 0;
    fh3MSubPtTrueDR[i]      = 0;
    fh3MTruePtTrueDR[i]     = 0;
    fh3PtTrueDeltaMDR[i]    = 0;
    fh3PtTrueDeltaMRelDR[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),
  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),
  fh3MSubPtTrueDR(0x0),
  fh3MTruePtTrueDR(0x0),
  fh3PtTrueDeltaMDR(0x0),
  fh3PtTrueDeltaMRelDR(0x0),
  fhnMassResponse(0x0)
{
  // Standard constructor.

  fh2MSubMatch         = new TH2F*[fNcentBins];
  fh2MSubPtRawAll      = new TH2F*[fNcentBins];
  fh3MSubPtRawDRMatch  = new TH3F*[fNcentBins];
  fh3MSubPtTrueDR      = new TH3F*[fNcentBins];
  fh3MTruePtTrueDR     = new TH3F*[fNcentBins];
  fh3PtTrueDeltaMDR    = new TH3F*[fNcentBins];
  fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
  fhnMassResponse      = new THnSparse*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh2MSubMatch[i]         = 0;
    fh2MSubPtRawAll[i]      = 0;
    fh3MSubPtRawDRMatch[i]  = 0;
    fh3MSubPtTrueDR[i]      = 0;
    fh3MTruePtTrueDR[i]     = 0;
    fh3PtTrueDeltaMDR[i]    = 0;
    fh3PtTrueDeltaMRelDR[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.;

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

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

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

    histName = Form("fh3PtTrueDeltaMDR_%d",i);
    histTitle = Form("fh3PtTrueDeltaMDR_%d;#it{p}_{T,true};#Delta %s",i,varName.Data());
    fh3PtTrueDeltaMDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDM,minDM,maxDM,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
    fOutput->Add(fh3PtTrueDeltaMDR[i]);

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

    histName = Form("fhnMassResponse_%d",i);
    histTitle = Form("fhnMassResponse_%d;%s sub;%s true;#it{p}_{T,sub};#it{p}_{T,true};#Delta R_{LJ}",i,varName.Data(),varName.Data());
    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);

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