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

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

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

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

#include "AliAODEvent.h"

#include "AliAnalysisTaskEmcalHJetMass.h"

ClassImp(EmcalHJetMassAnalysis::AliAnalysisTaskEmcalHJetMass)

namespace EmcalHJetMassAnalysis {

  //________________________________________________________________________
  AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass() : 
    AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalHJetMass", kTRUE),
    fContainerBase(0),
    fContainerUnsub(1),
    fMinFractionShared(0),
    fUseUnsubJet(0),
    fJetMassType(kRaw),
    fDPhiHJetMax(0.6),
    fTriggerTrackType(kInclusive),
    fPtTTMin(0),
    fPtTTMax(0),
    fRandom(0),
    fh1PtHadron(0),
    fh3PtHPtJDPhi(0),
    fh3PtJet1VsMassVsHPtAllSel(0),
    fh3PtJet1VsMassVsHPtAllSelMatch(0),
    fh3PtJet1VsMassVsHPtTagged(0),
    fh3PtJet1VsMassVsHPtTaggedMatch(0),
    fh3PtJet1VsRatVsHPtAllSel(0),
    fh3PtJet1VsRatVsHPtAllSelMatch(0),
    fh3PtJet1VsRatVsHPtTagged(0),
    fh3PtJet1VsRatVsHPtTaggedMatch(0)
  {
    // Default constructor.

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

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

    fPtTTMin = new TArrayF();
    fPtTTMax = new TArrayF();

    SetMakeGeneralHistograms(kTRUE);
  }

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

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

    fPtTTMin = new TArrayF();
    fPtTTMax = new TArrayF();

    SetMakeGeneralHistograms(kTRUE);
  }

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

    if(fRandom)      delete fRandom;
    if(fPtTTMin)     delete fPtTTMin;
    if(fPtTTMax)     delete fPtTTMax;
  }

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

    AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

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

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

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

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

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

    const Int_t nBinsPhi  = 18*4;
    const Double_t minPhi = -0.5*TMath::Pi();
    const Double_t maxPhi = 1.5*TMath::Pi();

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

      histName = TString::Format("fh3PtHPtJDPhi_%d",i);
      histTitle = TString::Format("%s;#it{p}_{T,h};#it{p}_{T,jet};#Delta#varphi_{h,jet}",histName.Data());
      fh3PtHPtJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPtH,minPtH,maxPtH,nBinsPt,minPt,maxPt,nBinsPhi,minPhi,maxPhi);
      fOutput->Add(fh3PtHPtJDPhi[i]);

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

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

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

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

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

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

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

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

    TH1::AddDirectory(oldStatus);

    fRandom = new TRandom3(0);

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

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

    AliParticleContainer *pCont = GetParticleContainer(0);
    AliJetContainer      *jCont = GetJetContainer(fContainerBase);
    if(!pCont || !jCont) return kFALSE;

    AliVParticle *vp = NULL;

    if(fTriggerTrackType==kInclusive) {
      pCont->ResetCurrentID();
      while((vp = pCont->GetNextAcceptParticle())) {
        fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all hadrons
        AliEmcalJet* jet = NULL;
        if(jCont) {
          jCont->ResetCurrentID();
          while((jet = jCont->GetNextAcceptJet())) {
            Double_t dphi = GetDeltaPhi(vp,jet)-TMath::Pi();
            fh3PtHPtJDPhi[fCentBin]->Fill(vp->Pt(),jet->Pt() - GetRhoVal(fContainerBase)*jet->Area(),dphi);
            if(dphi>fDPhiHJetMax) continue;
            FillHJetHistograms(vp->Pt(),jet);
          }
        }
      }
    }
    else if(fTriggerTrackType==kSingleInclusive) {
     for(Int_t it = 0; it<fPtTTMin->GetSize(); it++) {
       vp = GetSingleInclusiveTT(pCont,fPtTTMin->At(it),fPtTTMax->At(it));
       if(!vp) continue;
       fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all trigger tracks
       AliEmcalJet* jet = NULL;
       jCont->ResetCurrentID();
       while((jet = jCont->GetNextAcceptJet())) {
         Double_t dphi = GetDeltaPhi(vp,jet)-TMath::Pi();
         fh3PtHPtJDPhi[fCentBin]->Fill(vp->Pt(),jet->Pt() - GetRhoVal(fContainerBase)*jet->Area(),dphi);
         if(dphi>fDPhiHJetMax) continue;
         FillHJetHistograms(vp->Pt(),jet);
       }
     }//trigger track types
    }
    return kTRUE;
  }

  //________________________________________________________________________
  AliVParticle* AliAnalysisTaskEmcalHJetMass::GetSingleInclusiveTT(AliParticleContainer *pCont, Double_t ptmin, Double_t ptmax) const {
    AliVParticle *vp;
    TArrayI arr; arr.Set(pCont->GetNParticles());
    arr.Reset();
    Int_t counter = -1;
    pCont->ResetCurrentID();
    while((vp = pCont->GetNextAcceptParticle())) {
      if(vp->Pt()>=ptmin && vp->Pt()<ptmax ) {
        counter++;
        arr.SetAt(pCont->GetCurrentID(),counter);
      }
    }
    if(counter<0) return NULL;
    //select trigger track randomly
    fRandom->SetSeed(arr.At(0)); //random selection reproducible
    Double_t rnd = fRandom->Uniform() * counter;
    Int_t trigID = arr.At(TMath::FloorNint(rnd));
    vp = pCont->GetParticle(trigID);
    return vp;
  }

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

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

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

    if(jet->GetTagStatus()<1 || !jet->GetTaggedJet())
      return kFALSE;

    fh3PtJet1VsMassVsHPtTagged[fCentBin]->Fill(ptJet,mJet,pt);
    fh3PtJet1VsRatVsHPtTagged[fCentBin]->Fill(ptJet,rat,pt);

    if(fMinFractionShared>0. && fraction>fMinFractionShared) {
      fh3PtJet1VsMassVsHPtTaggedMatch[fCentBin]->Fill(ptJet,mJet,pt);
      fh3PtJet1VsRatVsHPtTaggedMatch[fCentBin]->Fill(ptJet,rat,pt);
    }
    return kTRUE;
  }

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

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

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


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

    return kTRUE;
  }

  //________________________________________________________________________
  void AliAnalysisTaskEmcalHJetMass::AddTriggerTrackPtCuts(Float_t min, Float_t max) {
    if(!fPtTTMin) fPtTTMin = new TArrayF();
    if(!fPtTTMax) fPtTTMax = new TArrayF();
    Int_t newSize = fPtTTMin->GetSize()+1;
    fPtTTMin->Set(newSize);
    fPtTTMax->Set(newSize);
    fPtTTMin->AddAt(min,newSize-1);
    fPtTTMax->AddAt(max,newSize-1);
  }

  //_______________________________________________________________________
  void AliAnalysisTaskEmcalHJetMass::Terminate(Option_t *) 
  {
    // Called once at the end of the analysis.
  }
}
Commit	Line	Data
1705bab1	1	//
	2	// Jet mass analysis task for jets recoiling from high pT hadron
	3	//
	4	// Author: M.Verweij
	5
	6	#include <TClonesArray.h>
	7	#include <TH1F.h>
	8	#include <TH2F.h>
	9	#include <TH3F.h>
	10	#include <THnSparse.h>
	11	#include <TList.h>
	12	#include <TLorentzVector.h>
	13	#include <TProfile.h>
	14	#include <TChain.h>
	15	#include <TSystem.h>
	16	#include <TFile.h>
	17	#include <TKey.h>
f887e746	18	#include <TArrayI.h>
	19	#include <TArrayF.h>
	20	#include <TRandom3.h>
1705bab1	21
	22	#include "AliVCluster.h"
	23	#include "AliVTrack.h"
	24	#include "AliEmcalJet.h"
	25	#include "AliRhoParameter.h"
	26	#include "AliLog.h"
	27	#include "AliEmcalParticle.h"
	28	#include "AliAnalysisManager.h"
	29	#include "AliJetContainer.h"
	30	#include "AliParticleContainer.h"
	31
	32	#include "AliAODEvent.h"
	33
	34	#include "AliAnalysisTaskEmcalHJetMass.h"
	35
f887e746	36	ClassImp(EmcalHJetMassAnalysis::AliAnalysisTaskEmcalHJetMass)
	37
	38	namespace EmcalHJetMassAnalysis {
	39
	40	//________________________________________________________________________
	41	AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass() :
	42	AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalHJetMass", kTRUE),
	43	fContainerBase(0),
	44	fContainerUnsub(1),
	45	fMinFractionShared(0),
	46	fUseUnsubJet(0),
	47	fJetMassType(kRaw),
	48	fDPhiHJetMax(0.6),
	49	fTriggerTrackType(kInclusive),
	50	fPtTTMin(0),
	51	fPtTTMax(0),
	52	fRandom(0),
	53	fh1PtHadron(0),
	54	fh3PtHPtJDPhi(0),
	55	fh3PtJet1VsMassVsHPtAllSel(0),
700200be	56	fh3PtJet1VsMassVsHPtAllSelMatch(0),
f887e746	57	fh3PtJet1VsMassVsHPtTagged(0),
	58	fh3PtJet1VsMassVsHPtTaggedMatch(0),
	59	fh3PtJet1VsRatVsHPtAllSel(0),
700200be	60	fh3PtJet1VsRatVsHPtAllSelMatch(0),
f887e746	61	fh3PtJet1VsRatVsHPtTagged(0),
	62	fh3PtJet1VsRatVsHPtTaggedMatch(0)
	63	{
	64	// Default constructor.
	65
	66	fh1PtHadron = new TH1F*[fNcentBins];
	67	fh3PtHPtJDPhi = new TH3F*[fNcentBins];
	68	fh3PtJet1VsMassVsHPtAllSel = new TH3F*[fNcentBins];
700200be	69	fh3PtJet1VsMassVsHPtAllSelMatch = new TH3F*[fNcentBins];
f887e746	70	fh3PtJet1VsMassVsHPtTagged = new TH3F*[fNcentBins];
	71	fh3PtJet1VsMassVsHPtTaggedMatch = new TH3F*[fNcentBins];
	72	fh3PtJet1VsRatVsHPtAllSel = new TH3F*[fNcentBins];
700200be	73	fh3PtJet1VsRatVsHPtAllSelMatch = new TH3F*[fNcentBins];
f887e746	74	fh3PtJet1VsRatVsHPtTagged = new TH3F*[fNcentBins];
	75	fh3PtJet1VsRatVsHPtTaggedMatch = new TH3F*[fNcentBins];
	76
	77	for (Int_t i = 0; i < fNcentBins; i++) {
	78	fh1PtHadron[i] = 0;
	79	fh3PtHPtJDPhi[i] = 0;
	80	fh3PtJet1VsMassVsHPtAllSel[i] = 0;
700200be	81	fh3PtJet1VsMassVsHPtAllSelMatch[i] = 0;
f887e746	82	fh3PtJet1VsMassVsHPtTagged[i] = 0;
	83	fh3PtJet1VsMassVsHPtTaggedMatch[i] = 0;
	84	fh3PtJet1VsRatVsHPtAllSel[i] = 0;
700200be	85	fh3PtJet1VsRatVsHPtAllSelMatch[i] = 0;
f887e746	86	fh3PtJet1VsRatVsHPtTagged[i] = 0;
	87	fh3PtJet1VsRatVsHPtTaggedMatch[i] = 0;
	88	}
1705bab1	89
f887e746	90	fPtTTMin = new TArrayF();
f887e746	91	fPtTTMax = new TArrayF();
1705bab1	92
f887e746	93	SetMakeGeneralHistograms(kTRUE);
1705bab1	94	}
1705bab1	95
f887e746	96	//________________________________________________________________________
	97	AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass(const char *name) :
	98	AliAnalysisTaskEmcalJet(name, kTRUE),
	99	fContainerBase(0),
	100	fContainerUnsub(1),
	101	fMinFractionShared(0),
	102	fUseUnsubJet(0),
	103	fJetMassType(kRaw),
	104	fDPhiHJetMax(0.6),
	105	fTriggerTrackType(kInclusive),
	106	fPtTTMin(0),
	107	fPtTTMax(0),
	108	fRandom(0),
	109	fh1PtHadron(0),
	110	fh3PtHPtJDPhi(0),
	111	fh3PtJet1VsMassVsHPtAllSel(0),
700200be	112	fh3PtJet1VsMassVsHPtAllSelMatch(0),
f887e746	113	fh3PtJet1VsMassVsHPtTagged(0),
	114	fh3PtJet1VsMassVsHPtTaggedMatch(0),
	115	fh3PtJet1VsRatVsHPtAllSel(0),
700200be	116	fh3PtJet1VsRatVsHPtAllSelMatch(0),
f887e746	117	fh3PtJet1VsRatVsHPtTagged(0),
	118	fh3PtJet1VsRatVsHPtTaggedMatch(0)
	119	{
	120	// Standard constructor.
	121
	122	fh1PtHadron = new TH1F*[fNcentBins];
	123	fh3PtHPtJDPhi = new TH3F*[fNcentBins];
	124	fh3PtJet1VsMassVsHPtAllSel = new TH3F*[fNcentBins];
700200be	125	fh3PtJet1VsMassVsHPtAllSelMatch = new TH3F*[fNcentBins];
f887e746	126	fh3PtJet1VsMassVsHPtTagged = new TH3F*[fNcentBins];
	127	fh3PtJet1VsMassVsHPtTaggedMatch = new TH3F*[fNcentBins];
	128	fh3PtJet1VsRatVsHPtAllSel = new TH3F*[fNcentBins];
700200be	129	fh3PtJet1VsRatVsHPtAllSelMatch = new TH3F*[fNcentBins];
f887e746	130	fh3PtJet1VsRatVsHPtTagged = new TH3F*[fNcentBins];
	131	fh3PtJet1VsRatVsHPtTaggedMatch = new TH3F*[fNcentBins];
	132
	133	for (Int_t i = 0; i < fNcentBins; i++) {
	134	fh1PtHadron[i] = 0;
	135	fh3PtHPtJDPhi[i] = 0;
	136	fh3PtJet1VsMassVsHPtAllSel[i] = 0;
700200be	137	fh3PtJet1VsMassVsHPtAllSelMatch[i] = 0;
f887e746	138	fh3PtJet1VsMassVsHPtTagged[i] = 0;
	139	fh3PtJet1VsMassVsHPtTaggedMatch[i] = 0;
	140	fh3PtJet1VsRatVsHPtAllSel[i] = 0;
700200be	141	fh3PtJet1VsRatVsHPtAllSelMatch[i] = 0;
f887e746	142	fh3PtJet1VsRatVsHPtTagged[i] = 0;
	143	fh3PtJet1VsRatVsHPtTaggedMatch[i] = 0;
	144	}
1705bab1	145
f887e746	146	fPtTTMin = new TArrayF();
f887e746	147	fPtTTMax = new TArrayF();
1705bab1	148
f887e746	149	SetMakeGeneralHistograms(kTRUE);
f887e746	150	}
ad3928e3	151
f887e746	152	//________________________________________________________________________
	153	AliAnalysisTaskEmcalHJetMass::~AliAnalysisTaskEmcalHJetMass()
	154	{
	155	// Destructor.
1705bab1	156
f887e746	157	if(fRandom) delete fRandom;
	158	if(fPtTTMin) delete fPtTTMin;
	159	if(fPtTTMax) delete fPtTTMax;
	160	}
ad3928e3	161
f887e746	162	//________________________________________________________________________
	163	void AliAnalysisTaskEmcalHJetMass::UserCreateOutputObjects()
	164	{
	165	// Create user output.
	166
	167	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
	168
	169	Bool_t oldStatus = TH1::AddDirectoryStatus();
	170	TH1::AddDirectory(kFALSE);
	171
	172	const Int_t nBinsPt = 200;
	173	const Double_t minPt = -50.;
	174	const Double_t maxPt = 150.;
	175
	176	const Int_t nBinsM = 100;
	177	const Double_t minM = -20.;
	178	const Double_t maxM = 80.;
	179
	180	const Int_t nBinsR = 100;
	181	const Double_t minR = -0.2;
	182	const Double_t maxR = 0.8;
	183
	184	const Int_t nBinsPtH = 100;
	185	const Double_t minPtH = 0.;
	186	const Double_t maxPtH = 100.;
	187
	188	const Int_t nBinsPhi = 18*4;
	189	const Double_t minPhi = -0.5*TMath::Pi();
	190	const Double_t maxPhi = 1.5*TMath::Pi();
	191
	192	TString histName = "";
	193	TString histTitle = "";
	194	for (Int_t i = 0; i < fNcentBins; i++) {
	195	histName = TString::Format("fh1PtHadron_%d",i);
	196	histTitle = TString::Format("%s;#it{p}_{T,h}",histName.Data());
	197	fh1PtHadron[i] = new TH1F(histName.Data(),histTitle.Data(),200.,0.,200.);
	198	fOutput->Add(fh1PtHadron[i]);
	199
	200	histName = TString::Format("fh3PtHPtJDPhi_%d",i);
	201	histTitle = TString::Format("%s;#it{p}_{T,h};#it{p}_{T,jet};#Delta#varphi_{h,jet}",histName.Data());
	202	fh3PtHPtJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPtH,minPtH,maxPtH,nBinsPt,minPt,maxPt,nBinsPhi,minPhi,maxPhi);
	203	fOutput->Add(fh3PtHPtJDPhi[i]);
	204
	205	histName = TString::Format("fh3PtJet1VsMassVsHPtAllSel_%d",i);
	206	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
	207	fh3PtJet1VsMassVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
	208	fOutput->Add(fh3PtJet1VsMassVsHPtAllSel[i]);
	209
700200be	210	histName = TString::Format("fh3PtJet1VsMassVsHPtAllSelMatch_%d",i);
	211	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
	212	fh3PtJet1VsMassVsHPtAllSelMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
	213	fOutput->Add(fh3PtJet1VsMassVsHPtAllSelMatch[i]);
	214
f887e746	215	histName = TString::Format("fh3PtJet1VsMassVsHPtTagged_%d",i);
	216	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
	217	fh3PtJet1VsMassVsHPtTagged[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
	218	fOutput->Add(fh3PtJet1VsMassVsHPtTagged[i]);
	219
	220	histName = TString::Format("fh3PtJet1VsMassVsHPtTaggedMatch_%d",i);
	221	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
	222	fh3PtJet1VsMassVsHPtTaggedMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
	223	fOutput->Add(fh3PtJet1VsMassVsHPtTaggedMatch[i]);
	224
	225	//
	226	histName = TString::Format("fh3PtJet1VsRatVsHPtAllSel_%d",i);
	227	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
	228	fh3PtJet1VsRatVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
	229	fOutput->Add(fh3PtJet1VsRatVsHPtAllSel[i]);
	230
700200be	231	histName = TString::Format("fh3PtJet1VsRatVsHPtAllSelMatch_%d",i);
	232	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
	233	fh3PtJet1VsRatVsHPtAllSelMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
	234	fOutput->Add(fh3PtJet1VsRatVsHPtAllSelMatch[i]);
	235
f887e746	236	histName = TString::Format("fh3PtJet1VsRatVsHPtTagged_%d",i);
	237	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
	238	fh3PtJet1VsRatVsHPtTagged[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
	239	fOutput->Add(fh3PtJet1VsRatVsHPtTagged[i]);
	240
	241	histName = TString::Format("fh3PtJet1VsRatVsHPtTaggedMatch_%d",i);
	242	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
	243	fh3PtJet1VsRatVsHPtTaggedMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
	244	fOutput->Add(fh3PtJet1VsRatVsHPtTaggedMatch[i]);
	245	}
1705bab1	246
f887e746	247	TH1::AddDirectory(oldStatus);
1705bab1	248
f887e746	249	fRandom = new TRandom3(0);
1705bab1	250
f887e746	251	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
1705bab1	252	}
1705bab1	253
f887e746	254	//________________________________________________________________________
	255	Bool_t AliAnalysisTaskEmcalHJetMass::Run()
	256	{
	257	// Run analysis code here, if needed. It will be executed before FillHistograms().
	258
	259	AliParticleContainer *pCont = GetParticleContainer(0);
	260	AliJetContainer *jCont = GetJetContainer(fContainerBase);
	261	if(!pCont \|\| !jCont) return kFALSE;
	262
	263	AliVParticle *vp = NULL;
	264
	265	if(fTriggerTrackType==kInclusive) {
	266	pCont->ResetCurrentID();
	267	while((vp = pCont->GetNextAcceptParticle())) {
	268	fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all hadrons
	269	AliEmcalJet* jet = NULL;
	270	if(jCont) {
	271	jCont->ResetCurrentID();
	272	while((jet = jCont->GetNextAcceptJet())) {
	273	Double_t dphi = GetDeltaPhi(vp,jet)-TMath::Pi();
	274	fh3PtHPtJDPhi[fCentBin]->Fill(vp->Pt(),jet->Pt() - GetRhoVal(fContainerBase)*jet->Area(),dphi);
	275	if(dphi>fDPhiHJetMax) continue;
	276	FillHJetHistograms(vp->Pt(),jet);
	277	}
	278	}
	279	}
	280	}
	281	else if(fTriggerTrackType==kSingleInclusive) {
700200be	282	for(Int_t it = 0; it<fPtTTMin->GetSize(); it++) {
	283	vp = GetSingleInclusiveTT(pCont,fPtTTMin->At(it),fPtTTMax->At(it));
	284	if(!vp) continue;
	285	fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all trigger tracks
	286	AliEmcalJet* jet = NULL;
	287	jCont->ResetCurrentID();
	288	while((jet = jCont->GetNextAcceptJet())) {
	289	Double_t dphi = GetDeltaPhi(vp,jet)-TMath::Pi();
	290	fh3PtHPtJDPhi[fCentBin]->Fill(vp->Pt(),jet->Pt() - GetRhoVal(fContainerBase)*jet->Area(),dphi);
	291	if(dphi>fDPhiHJetMax) continue;
	292	FillHJetHistograms(vp->Pt(),jet);
	293	}
	294	}//trigger track types
1705bab1	295	}
f887e746	296	return kTRUE;
1705bab1	297	}
1705bab1	298
700200be	299	//________________________________________________________________________
	300	AliVParticle* AliAnalysisTaskEmcalHJetMass::GetSingleInclusiveTT(AliParticleContainer *pCont, Double_t ptmin, Double_t ptmax) const {
	301	AliVParticle *vp;
	302	TArrayI arr; arr.Set(pCont->GetNParticles());
	303	arr.Reset();
	304	Int_t counter = -1;
	305	pCont->ResetCurrentID();
	306	while((vp = pCont->GetNextAcceptParticle())) {
	307	if(vp->Pt()>=ptmin && vp->Pt()<ptmax ) {
	308	counter++;
	309	arr.SetAt(pCont->GetCurrentID(),counter);
	310	}
	311	}
	312	if(counter<0) return NULL;
	313	//select trigger track randomly
	314	fRandom->SetSeed(arr.At(0)); //random selection reproducible
	315	Double_t rnd = fRandom->Uniform() * counter;
	316	Int_t trigID = arr.At(TMath::FloorNint(rnd));
	317	vp = pCont->GetParticle(trigID);
	318	return vp;
	319	}
	320
f887e746	321	//________________________________________________________________________
	322	Bool_t AliAnalysisTaskEmcalHJetMass::FillHJetHistograms(Double_t pt, const AliEmcalJet *jet)
	323	{
	324	// Fill hadron-jet histograms.
	325	Double_t ptJet = jet->Pt() - GetRhoVal(fContainerBase)*jet->Area();
	326	Double_t mJet = GetJetMass(jet);
	327	Double_t rat = -1.;
	328	if(ptJet<0. \|\| ptJet>0.) rat = mJet/ptJet;
	329
	330	fh3PtJet1VsMassVsHPtAllSel[fCentBin]->Fill(ptJet,mJet,pt);
	331	fh3PtJet1VsRatVsHPtAllSel[fCentBin]->Fill(ptJet,rat,pt);
	332
f887e746	333	Double_t fraction = -1.;
	334	if(fUseUnsubJet) {
	335	AliEmcalJet *jetUS = NULL;
	336	AliJetContainer *jetContUS = GetJetContainer(fContainerUnsub);
	337	Int_t ifound = 0;
	338	Int_t ilab = -1;
	339	for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
	340	jetUS = jetContUS->GetJet(i);
	341	if(jetUS->GetLabel()==jet->GetLabel()) {
	342	ifound++;
	343	if(ifound==1) ilab = i;
	344	}
1705bab1	345	}
f887e746	346	if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
	347	if(ifound==0) jetUS = 0x0;
	348	else jetUS = jetContUS->GetJet(ilab);
	349	fraction = jetContUS->GetFractionSharedPt(jetUS);
	350	} else {
	351	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
	352	fraction = jetCont->GetFractionSharedPt(jet);
1705bab1	353	}
1705bab1	354
700200be	355	if(fMinFractionShared>0. && fraction>fMinFractionShared) {
	356	fh3PtJet1VsMassVsHPtAllSelMatch[fCentBin]->Fill(ptJet,mJet,pt);
	357	fh3PtJet1VsRatVsHPtAllSelMatch[fCentBin]->Fill(ptJet,rat,pt);
	358	}
	359
	360	if(jet->GetTagStatus()<1 \|\| !jet->GetTaggedJet())
	361	return kFALSE;
	362
	363	fh3PtJet1VsMassVsHPtTagged[fCentBin]->Fill(ptJet,mJet,pt);
	364	fh3PtJet1VsRatVsHPtTagged[fCentBin]->Fill(ptJet,rat,pt);
	365
f887e746	366	if(fMinFractionShared>0. && fraction>fMinFractionShared) {
	367	fh3PtJet1VsMassVsHPtTaggedMatch[fCentBin]->Fill(ptJet,mJet,pt);
	368	fh3PtJet1VsRatVsHPtTaggedMatch[fCentBin]->Fill(ptJet,rat,pt);
	369	}
	370	return kTRUE;
1705bab1	371	}
1705bab1	372
f887e746	373	//________________________________________________________________________
	374	Double_t AliAnalysisTaskEmcalHJetMass::GetJetMass(const AliEmcalJet *jet) const {
	375	//calc subtracted jet mass
	376	if(fJetMassType==kRaw)
	377	return jet->M();
	378	else if(fJetMassType==kDeriv)
	379	return jet->GetSecondOrderSubtracted();
1705bab1	380
f887e746	381	return 0;
f887e746	382	}
1705bab1	383
f887e746	384	//________________________________________________________________________
	385	Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(const AliVParticle vp, const AliEmcalJet jet) const {
	386	// Calculate azimuthal angle between particle and jet. range:[-0.5\pi,1.5\pi]
	387	return GetDeltaPhi(vp->Phi(),jet->Phi());
	388	}
1705bab1	389
f887e746	390	//________________________________________________________________________
	391	Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(Double_t phi1, Double_t phi2) const {
	392	// Calculate azimuthal angle between phi1 and phi2. range:[-0.5\pi,1.5\pi]
	393	Double_t dPhi = phi1-phi2;
	394	if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
	395	if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
	396	return dPhi;
	397	}
1705bab1	398
1705bab1	399
f887e746	400	//________________________________________________________________________
	401	Bool_t AliAnalysisTaskEmcalHJetMass::RetrieveEventObjects() {
	402	//
	403	// retrieve event objects
	404	//
	405	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
	406	return kFALSE;
	407
	408	return kTRUE;
	409	}
1705bab1	410
f887e746	411	//________________________________________________________________________
	412	void AliAnalysisTaskEmcalHJetMass::AddTriggerTrackPtCuts(Float_t min, Float_t max) {
	413	if(!fPtTTMin) fPtTTMin = new TArrayF();
	414	if(!fPtTTMax) fPtTTMax = new TArrayF();
	415	Int_t newSize = fPtTTMin->GetSize()+1;
	416	fPtTTMin->Set(newSize);
	417	fPtTTMax->Set(newSize);
	418	fPtTTMin->AddAt(min,newSize-1);
	419	fPtTTMax->AddAt(max,newSize-1);
	420	}
1705bab1	421
f887e746	422	//_______________________________________________________________________
	423	void AliAnalysisTaskEmcalHJetMass::Terminate(Option_t *)
	424	{
	425	// Called once at the end of the analysis.
	426	}
1705bab1	427	}
1705bab1	428