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

//
// Jet tagger analysis task.
//
// Author: M.Verweij

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

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

#include "AliAODEvent.h"
#include "AliESDEvent.h"

#include "AliAnalysisTaskEmcalJetTagger.h"

ClassImp(AliAnalysisTaskEmcalJetTagger)

//________________________________________________________________________
AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTagger", kTRUE),
  fJetTaggingType(kTag),
  fJetTaggingMethod(kGeo),
  fContainerBase(0),
  fContainerTag(1),
  fMinFractionShared(0),
  fMatchingDone(0),
  fh3PtJet1VsDeltaEtaDeltaPhi(0),
  fh2PtJet1VsDeltaR(0),
  fh2PtJet2VsFraction(0),
  fh2PtJet1VsLeadPtAllSel(0),
  fh2PtJet1VsLeadPtTagged(0),
  fh2PtJet1VsPtJet2(0),
  fh3PtJetDEtaDPhiConst(0),
  fh2PtJetDRConst(0),
  fh3PtJetAreaDRConst(0)
{
  // Default constructor.

  fh3PtJet1VsDeltaEtaDeltaPhi  = new TH3F*[fNcentBins];
  fh2PtJet1VsDeltaR            = new TH2F*[fNcentBins];
  fh2PtJet2VsFraction          = new TH2F*[fNcentBins];
  fh2PtJet1VsLeadPtAllSel      = new TH2F*[fNcentBins];
  fh2PtJet1VsLeadPtTagged      = new TH2F*[fNcentBins];
  fh2PtJet1VsPtJet2            = new TH2F*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
    fh2PtJet1VsDeltaR[i]           = 0;
    fh2PtJet2VsFraction[i]         = 0;
    fh2PtJet1VsLeadPtAllSel[i]     = 0;
    fh2PtJet1VsLeadPtTagged[i]     = 0;
    fh2PtJet1VsPtJet2[i]           = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
  
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),  
  fJetTaggingType(kTag),
  fJetTaggingMethod(kGeo),
  fContainerBase(0),
  fContainerTag(1),
  fMinFractionShared(0),
  fMatchingDone(0),
  fh3PtJet1VsDeltaEtaDeltaPhi(0),
  fh2PtJet1VsDeltaR(0),
  fh2PtJet2VsFraction(0),
  fh2PtJet1VsLeadPtAllSel(0),
  fh2PtJet1VsLeadPtTagged(0),
  fh2PtJet1VsPtJet2(0),
  fh3PtJetDEtaDPhiConst(0),
  fh2PtJetDRConst(0),
  fh3PtJetAreaDRConst(0)
{
  // Standard constructor.

  fh3PtJet1VsDeltaEtaDeltaPhi = new TH3F*[fNcentBins];
  fh2PtJet1VsDeltaR           = new TH2F*[fNcentBins];
  fh2PtJet2VsFraction         = new TH2F*[fNcentBins];
  fh2PtJet1VsLeadPtAllSel     = new TH2F*[fNcentBins];
  fh2PtJet1VsLeadPtTagged     = new TH2F*[fNcentBins];
  fh2PtJet1VsPtJet2           = new TH2F*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
    fh2PtJet1VsDeltaR[i]           = 0;
    fh2PtJet2VsFraction[i]         = 0;
    fh2PtJet1VsLeadPtAllSel[i]     = 0;
    fh2PtJet1VsLeadPtTagged[i]     = 0;
    fh2PtJet1VsPtJet2[i]           = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

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

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

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

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

  const Int_t nBinsPt          = 250;
  const Int_t nBinsDPhi        = 100;
  const Int_t nBinsDEta        = 100;
  const Int_t nBinsDR          = 50;
  const Int_t nBinsFraction    = 101;

  const Double_t minPt       = -50.;
  const Double_t maxPt       = 200.;
  const Double_t minDPhi     = -0.5;
  const Double_t maxDPhi     =  0.5;
  const Double_t minDEta     = -0.5;
  const Double_t maxDEta     =  0.5;
  const Double_t minDR       =  0.;
  const Double_t maxDR       =  0.5;
  const Double_t minFraction =  -0.005;
  const Double_t maxFraction =  1.005;

  TString histName = "";
  TString histTitle = "";
  for (Int_t i = 0; i < fNcentBins; i++) {

    histName = TString::Format("fh3PtJet1VsDeltaEtaDeltaPhi_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta#eta};#it{#Delta#varphi}",histName.Data());
    fh3PtJet1VsDeltaEtaDeltaPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDEta,minDEta,maxDEta,nBinsDPhi,minDPhi,maxDPhi);
    fOutput->Add(fh3PtJet1VsDeltaEtaDeltaPhi[i]);

    histName = TString::Format("fh2PtJet1VsDeltaR_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta R}",histName.Data());
    fh2PtJet1VsDeltaR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDR,minDR,maxDR);
    fOutput->Add(fh2PtJet1VsDeltaR[i]);

    histName = TString::Format("fh2PtJet2VsFraction_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet2};#it{f}_{shared}",histName.Data());
    fh2PtJet2VsFraction[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsFraction,minFraction,maxFraction);
    fOutput->Add(fh2PtJet2VsFraction[i]);

    histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
    fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
    fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);

    histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
    fh2PtJet1VsLeadPtTagged[i] =  new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
    fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);

    histName = TString::Format("fh2PtJet1VsPtJet2_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,jet2}",histName.Data());
    fh2PtJet1VsPtJet2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
    fOutput->Add(fh2PtJet1VsPtJet2[i]);

  }

  fh3PtJetDEtaDPhiConst = new TH3F("fh3PtJetDEtaDPhiConst","fh3PtJetDEtaDPhiConst;pT;#Delta #eta;#Delta #varphi",nBinsPt,minPt,maxPt,nBinsDEta,-1.,1.,nBinsDPhi,-1.,1.);
  fOutput->Add(fh3PtJetDEtaDPhiConst);

  fh2PtJetDRConst = new TH2F("fh2PtJetDRConst","fh2PtJetDRConst;pT;#Delta R",nBinsPt,minPt,maxPt,100,0.,1.);
  fOutput->Add(fh2PtJetDRConst);

  fh3PtJetAreaDRConst = new TH3F("fh3PtJetAreaDRConst","fh3PtJetAreaDRConst;pT;A;#Delta R",nBinsPt,minPt,maxPt,100,0.,1.,100,0.,1.);
  fOutput->Add(fh3PtJetAreaDRConst);

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

  TH1::AddDirectory(oldStatus);

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

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

  MatchJetsGeo(fContainerBase,fContainerTag,0,0.3,2);

  //  if(fJetTaggingMethod==kFraction)

  return kTRUE;
}

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

  AliEmcalJet *jet1 = NULL;
  AliJetContainer *jetCont = GetJetContainer(fContainerBase);
  if(!jetCont) return kFALSE;
  jetCont->ResetCurrentID();
  while((jet1 = jetCont->GetNextAcceptJet())) {
    Double_t ptJet1 =  jet1->Pt() - jetCont->GetRhoVal()*jet1->Area();
    fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());

    //fill histo with angle between jet axis and constituents
    for(Int_t icc=0; icc<jet1->GetNumberOfTracks(); icc++) {
      AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(icc, jetCont->GetParticleContainer()->GetArray()));//fTracks));
      if(!vp) continue;
      Double_t dEta = jet1->Eta()-vp->Eta();
      Double_t dPhi = jet1->Phi()-vp->Phi();
      if(dPhi<TMath::Pi()) dPhi+=TMath::TwoPi();
      if(dPhi>TMath::Pi()) dPhi-=TMath::TwoPi();
      fh3PtJetDEtaDPhiConst->Fill(ptJet1,dEta,dPhi);

      Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
      fh2PtJetDRConst->Fill(ptJet1,dR);
      fh3PtJetAreaDRConst->Fill(ptJet1,jet1->Area(),dR);
    }

    if(jet1->GetTagStatus()<1 && fJetTaggingType==kTag)
      continue;

    AliEmcalJet *jet2 = NULL;
    if(fJetTaggingType==kTag)     jet2 = jet1->GetTaggedJet();
    if(fJetTaggingType==kClosest) jet2 = jet1->ClosestJet();
    if(!jet2) continue;

    Double_t ptJet2 =  jet2->Pt() - GetRhoVal(fContainerTag)*jet2->Area();
    Double_t fraction = jetCont->GetFractionSharedPt(jet1);
    fh2PtJet2VsFraction[fCentBin]->Fill(ptJet2,fraction);

    if(fraction<fMinFractionShared && fJetTaggingType==kClosest)
      continue;

    fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
    fh2PtJet1VsPtJet2[fCentBin]->Fill(ptJet1,ptJet2);

    Double_t dPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
    if(dPhi>TMath::Pi())
      dPhi -= TMath::TwoPi();
    if(dPhi<(-1.*TMath::Pi()))
      dPhi += TMath::TwoPi();  
    
    fh3PtJet1VsDeltaEtaDeltaPhi[fCentBin]->Fill(ptJet1,jet1->Eta()-jet2->Eta(),dPhi);
    fh2PtJet1VsDeltaR[fCentBin]->Fill(ptJet1,GetDeltaR(jet1,jet2));
  }
  return kTRUE;
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetTagger::ResetTagging(const Int_t c) {

  //Reset tagging of container c

  for(int i = 0;i<GetNJets(c);i++){
    AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetJetFromArray(i, c));
    if(fJetTaggingType==kClosest)
      jet->ResetMatching();
    else if(fJetTaggingType==kTag) {
      jet->SetTaggedJet(0x0);
      jet->SetTagStatus(-1);
    }
  }
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
						 Int_t iDebug, Float_t maxDist, Int_t type, Bool_t bReset) {

  //
  // Match the full jets to the corresponding charged jets
  // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
  // type: 0 = use acceptance cuts of container  1 = allow 0.1 one more for c2 in eta 2 = allow 0.1 more in eta and phi for c2

  if(c1<0) c1 = fContainerBase;
  if(c2<0) c2 = fContainerTag;

  const Int_t nJets1 = GetNJets(c1);
  const Int_t nJets2 = GetNJets(c2);

  if(nJets1==0 || nJets2==0) return;

  if(bReset) {
    ResetTagging(c1);
    ResetTagging(c2);
  }
  fMatchingDone = kFALSE;

  TArrayI faMatchIndex1;
  faMatchIndex1.Set(nJets2+1);
  faMatchIndex1.Reset(-1);

  TArrayI faMatchIndex2;
  faMatchIndex2.Set(nJets1+1);
  faMatchIndex2.Reset(-1);

  static TArrayS iFlag(nJets1*nJets2);
  if(iFlag.GetSize()<(nJets1*nJets2)){
    iFlag.Set(nJets1*nJets2+1);
  }
  iFlag.Reset(0);


  AliJetContainer *cont2 = GetJetContainer(c2);

  // find the closest distance to the full jet
  for(int i = 0;i<nJets1;i++){

    AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
    if(!jet1) continue;

    Float_t dist = maxDist;
    
    for(int j = 0;j <nJets2; j++){
      AliEmcalJet *jet2 = 0x0;
      if(type==0)
	jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
      else {
	jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
	if(!jet2) continue;
	if(type>0) {
	  if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) || jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
	    continue;
	  if(type==2) {
	    if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) || jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
	      continue;
	  }
	}
      }
      if(!jet2)
	continue;

      Double_t dR = GetDeltaR(jet1,jet2);
      if(dR<dist && dR<maxDist){
	faMatchIndex2[i]=j;
	dist = dR;
      }
    }//j jet loop
    if(faMatchIndex2[i]>=0) {
      iFlag[i*nJets2+faMatchIndex2[i]]+=1;//j closest to i
      if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",i,faMatchIndex2[i],dist,i*nJets2+faMatchIndex2[i],iFlag[i*nJets2+faMatchIndex2[i]]);
    }
  
  }//i jet loop


  // other way around
  for(int j = 0;j<nJets2;j++){
    AliEmcalJet *jet2 = 0x0;
    if(type==0)
      jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
    else {
      jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
      if(!jet2) continue;;
      if(type>0) {
	if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) || jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
	  continue;
	if(type==2) {
	  if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) || jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
	    continue;
	}
      }
    }
    if(!jet2)
      continue;

    Float_t dist = maxDist;
    for(int i = 0;i<nJets1;i++){
      AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
      if(!jet1)
	continue;

      Double_t dR = GetDeltaR(jet1,jet2); 
      if(dR<dist && dR<maxDist){
	faMatchIndex1[j]=i;
        dist = dR;
      }   
    }
    if(faMatchIndex1[j]>=0) {
      iFlag[faMatchIndex1[j]*nJets2+j]+=2;//i closest to j
      if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faMatchIndex1[j],j,dist,faMatchIndex1[j]*nJets2+j,iFlag[faMatchIndex1[j]*nJets2+j]);
    }

  }
    
  // check for "true" correlations
  for(int i = 0;i<nJets1;i++){
    AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetJetFromArray(i, c1));
    for(int j = 0;j<nJets2;j++){
      AliEmcalJet *jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
      AliDebug(11,Form("%s: Flag[%d][%d] %d ",GetName(),i,j,iFlag[i*nJets2+j]));
      
      // we have a uniqe correlation
      if(iFlag[i*nJets2+j]==3){
	
	Double_t dR = GetDeltaR(jet1,jet2); 
	if(iDebug>1) Printf("closest jets %d  %d  dR =  %f",j,i,dR);
	
	if(fJetTaggingType==kTag) {
	  jet1->SetTaggedJet(jet2);
	  jet1->SetTagStatus(1);
	  
	  jet2->SetTaggedJet(jet1);
	  jet2->SetTagStatus(1);
	}
	else if(fJetTaggingType==kClosest) {
	  jet1->SetClosestJet(jet2,dR);
	  jet2->SetClosestJet(jet1,dR);
	}
      }
    }
  }
  fMatchingDone = kTRUE;
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
  //
  // Helper function to calculate the distance between two jets
  //

  Double_t dPhi = jet1->Phi() - jet2->Phi();
  if(dPhi>TMath::Pi())
    dPhi -= TMath::TwoPi();
  if(dPhi<(-1.*TMath::Pi()))
    dPhi += TMath::TwoPi();
  Double_t dEta = jet1->Eta() - jet2->Eta();
  Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
  return dR;
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
  //
  // Calculate azimuthal angle between the axises of the jets
  //

  return GetDeltaPhi(jet1->Phi(),jet2->Phi());

}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(Double_t phi1,Double_t phi2) {
  //
  // Calculate azimuthal angle between the axises of the jets
  //

  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;
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalJetTagger::GetFractionSharedPt(const AliEmcalJet *jet1, const AliEmcalJet *jet2) const
{
  //
  // Get fraction of shared pT between matched full and charged jet
  // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
  //

  Double_t fraction = 0.;
  Double_t jetPt2 = jet2->Pt();
 
  if(jetPt2>0) {

    AliDebug(2,Form("%s: nConstituents_jet1: %d, nConstituents_jet2: %d, tracks_jet1: %d  tracks_jet2: %d clusters_jet1: %d clusters_jet2: %d",GetName(),jet1->GetNumberOfConstituents(),jet2->GetNumberOfConstituents(),jet1->GetNumberOfTracks(),jet2->GetNumberOfTracks(),jet1->GetNumberOfClusters(),jet2->GetNumberOfClusters()));
    
    Double_t sumPt = 0.;
    AliVParticle *vpf = 0x0;
    Int_t iFound = 0;
    for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
      Int_t idx = (Int_t)jet2->TrackAt(icc);
      iFound = 0;
      for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
	if(idx == jet1->TrackAt(icf) && iFound==0 ) {
	  iFound=1;
	  vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fTracks));
	  if(vpf) sumPt += vpf->Pt();
	  continue;
	}
      }
    }
    fraction = sumPt/jetPt2;
  } else fraction = -1.;
  return fraction;
}

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

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

  return kTRUE;

}

//_______________________________________________________________________
void AliAnalysisTaskEmcalJetTagger::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
8612dfc8	1	//
	2	// Jet tagger analysis task.
	3	//
	4	// Author: M.Verweij
	5
	6	#include <TClonesArray.h>
	7	#include <TH1F.h>
	8	#include <TH2F.h>
	9	#include <TH3F.h>
	10	#include <THnSparse.h>
	11	#include <TList.h>
	12	#include <TLorentzVector.h>
	13	#include <TProfile.h>
	14	#include <TChain.h>
	15	#include <TSystem.h>
	16	#include <TFile.h>
	17	#include <TKey.h>
	18
	19	#include "AliVCluster.h"
	20	#include "AliVTrack.h"
	21	#include "AliEmcalJet.h"
	22	#include "AliRhoParameter.h"
	23	#include "AliLog.h"
	24	#include "AliEmcalParticle.h"
	25	#include "AliMCEvent.h"
	26	#include "AliGenPythiaEventHeader.h"
	27	#include "AliAODMCHeader.h"
	28	#include "AliMCEvent.h"
	29	#include "AliAnalysisManager.h"
	30	#include "AliJetContainer.h"
ba7663ad	31	#include "AliParticleContainer.h"
8612dfc8	32
396ffc82 MV	33	#include "AliAODEvent.h"
	34	#include "AliESDEvent.h"
	35
8612dfc8	36	#include "AliAnalysisTaskEmcalJetTagger.h"
	37
	38	ClassImp(AliAnalysisTaskEmcalJetTagger)
	39
	40	//________________________________________________________________________
	41	AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger() :
	42	AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTagger", kTRUE),
	43	fJetTaggingType(kTag),
	44	fJetTaggingMethod(kGeo),
	45	fContainerBase(0),
	46	fContainerTag(1),
ba7663ad	47	fMinFractionShared(0),
8612dfc8	48	fMatchingDone(0),
	49	fh3PtJet1VsDeltaEtaDeltaPhi(0),
	50	fh2PtJet1VsDeltaR(0),
ba7663ad	51	fh2PtJet2VsFraction(0),
8612dfc8	52	fh2PtJet1VsLeadPtAllSel(0),
8612dfc8	53	fh2PtJet1VsLeadPtTagged(0),
396ffc82 MV	54	fh2PtJet1VsPtJet2(0),
	55	fh3PtJetDEtaDPhiConst(0),
	56	fh2PtJetDRConst(0),
	57	fh3PtJetAreaDRConst(0)
8612dfc8	58	{
	59	// Default constructor.
	60
	61	fh3PtJet1VsDeltaEtaDeltaPhi = new TH3F*[fNcentBins];
	62	fh2PtJet1VsDeltaR = new TH2F*[fNcentBins];
ba7663ad	63	fh2PtJet2VsFraction = new TH2F*[fNcentBins];
8612dfc8	64	fh2PtJet1VsLeadPtAllSel = new TH2F*[fNcentBins];
	65	fh2PtJet1VsLeadPtTagged = new TH2F*[fNcentBins];
	66	fh2PtJet1VsPtJet2 = new TH2F*[fNcentBins];
	67
	68	for (Int_t i = 0; i < fNcentBins; i++) {
	69	fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
	70	fh2PtJet1VsDeltaR[i] = 0;
ba7663ad	71	fh2PtJet2VsFraction[i] = 0;
8612dfc8	72	fh2PtJet1VsLeadPtAllSel[i] = 0;
	73	fh2PtJet1VsLeadPtTagged[i] = 0;
	74	fh2PtJet1VsPtJet2[i] = 0;
	75	}
	76
	77	SetMakeGeneralHistograms(kTRUE);
	78
	79	}
	80
	81	//________________________________________________________________________
	82	AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger(const char *name) :
	83	AliAnalysisTaskEmcalJet(name, kTRUE),
	84	fJetTaggingType(kTag),
	85	fJetTaggingMethod(kGeo),
	86	fContainerBase(0),
	87	fContainerTag(1),
ba7663ad	88	fMinFractionShared(0),
8612dfc8	89	fMatchingDone(0),
	90	fh3PtJet1VsDeltaEtaDeltaPhi(0),
	91	fh2PtJet1VsDeltaR(0),
ba7663ad	92	fh2PtJet2VsFraction(0),
8612dfc8	93	fh2PtJet1VsLeadPtAllSel(0),
8612dfc8	94	fh2PtJet1VsLeadPtTagged(0),
396ffc82 MV	95	fh2PtJet1VsPtJet2(0),
	96	fh3PtJetDEtaDPhiConst(0),
	97	fh2PtJetDRConst(0),
	98	fh3PtJetAreaDRConst(0)
8612dfc8	99	{
	100	// Standard constructor.
	101
	102	fh3PtJet1VsDeltaEtaDeltaPhi = new TH3F*[fNcentBins];
	103	fh2PtJet1VsDeltaR = new TH2F*[fNcentBins];
ba7663ad	104	fh2PtJet2VsFraction = new TH2F*[fNcentBins];
8612dfc8	105	fh2PtJet1VsLeadPtAllSel = new TH2F*[fNcentBins];
	106	fh2PtJet1VsLeadPtTagged = new TH2F*[fNcentBins];
	107	fh2PtJet1VsPtJet2 = new TH2F*[fNcentBins];
	108
	109	for (Int_t i = 0; i < fNcentBins; i++) {
	110	fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
	111	fh2PtJet1VsDeltaR[i] = 0;
ba7663ad	112	fh2PtJet2VsFraction[i] = 0;
8612dfc8	113	fh2PtJet1VsLeadPtAllSel[i] = 0;
	114	fh2PtJet1VsLeadPtTagged[i] = 0;
	115	fh2PtJet1VsPtJet2[i] = 0;
	116	}
	117
	118	SetMakeGeneralHistograms(kTRUE);
	119	}
	120
	121	//________________________________________________________________________
	122	AliAnalysisTaskEmcalJetTagger::~AliAnalysisTaskEmcalJetTagger()
	123	{
	124	// Destructor.
	125	}
	126
	127	//________________________________________________________________________
	128	void AliAnalysisTaskEmcalJetTagger::UserCreateOutputObjects()
	129	{
	130	// Create user output.
	131
	132	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
	133
	134	Bool_t oldStatus = TH1::AddDirectoryStatus();
	135	TH1::AddDirectory(kFALSE);
	136
	137	const Int_t nBinsPt = 250;
	138	const Int_t nBinsDPhi = 100;
	139	const Int_t nBinsDEta = 100;
	140	const Int_t nBinsDR = 50;
ba7663ad	141	const Int_t nBinsFraction = 101;
	142
	143	const Double_t minPt = -50.;
	144	const Double_t maxPt = 200.;
	145	const Double_t minDPhi = -0.5;
	146	const Double_t maxDPhi = 0.5;
	147	const Double_t minDEta = -0.5;
	148	const Double_t maxDEta = 0.5;
	149	const Double_t minDR = 0.;
	150	const Double_t maxDR = 0.5;
	151	const Double_t minFraction = -0.005;
	152	const Double_t maxFraction = 1.005;
8612dfc8	153
	154	TString histName = "";
	155	TString histTitle = "";
	156	for (Int_t i = 0; i < fNcentBins; i++) {
	157
	158	histName = TString::Format("fh3PtJet1VsDeltaEtaDeltaPhi_%d",i);
	159	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta#eta};#it{#Delta#varphi}",histName.Data());
	160	fh3PtJet1VsDeltaEtaDeltaPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDEta,minDEta,maxDEta,nBinsDPhi,minDPhi,maxDPhi);
	161	fOutput->Add(fh3PtJet1VsDeltaEtaDeltaPhi[i]);
	162
	163	histName = TString::Format("fh2PtJet1VsDeltaR_%d",i);
	164	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta R}",histName.Data());
	165	fh2PtJet1VsDeltaR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDR,minDR,maxDR);
	166	fOutput->Add(fh2PtJet1VsDeltaR[i]);
	167
ba7663ad	168	histName = TString::Format("fh2PtJet2VsFraction_%d",i);
	169	histTitle = TString::Format("%s;#it{p}_{T,jet2};#it{f}_{shared}",histName.Data());
	170	fh2PtJet2VsFraction[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsFraction,minFraction,maxFraction);
	171	fOutput->Add(fh2PtJet2VsFraction[i]);
	172
8612dfc8	173	histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
	174	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
	175	fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
	176	fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
	177
	178	histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
	179	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
	180	fh2PtJet1VsLeadPtTagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
	181	fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
	182
	183	histName = TString::Format("fh2PtJet1VsPtJet2_%d",i);
	184	histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,jet2}",histName.Data());
	185	fh2PtJet1VsPtJet2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
	186	fOutput->Add(fh2PtJet1VsPtJet2[i]);
	187
	188	}
	189
396ffc82 MV	190	fh3PtJetDEtaDPhiConst = new TH3F("fh3PtJetDEtaDPhiConst","fh3PtJetDEtaDPhiConst;pT;#Delta #eta;#Delta #varphi",nBinsPt,minPt,maxPt,nBinsDEta,-1.,1.,nBinsDPhi,-1.,1.);
	191	fOutput->Add(fh3PtJetDEtaDPhiConst);
	192
	193	fh2PtJetDRConst = new TH2F("fh2PtJetDRConst","fh2PtJetDRConst;pT;#Delta R",nBinsPt,minPt,maxPt,100,0.,1.);
	194	fOutput->Add(fh2PtJetDRConst);
	195
	196	fh3PtJetAreaDRConst = new TH3F("fh3PtJetAreaDRConst","fh3PtJetAreaDRConst;pT;A;#Delta R",nBinsPt,minPt,maxPt,100,0.,1.,100,0.,1.);
	197	fOutput->Add(fh3PtJetAreaDRConst);
	198
8612dfc8	199	// =========== Switch on Sumw2 for all histos ===========
	200	for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
	201	TH1 h1 = dynamic_cast<TH1>(fOutput->At(i));
	202	if (h1){
	203	h1->Sumw2();
	204	continue;
	205	}
	206	THnSparse hn = dynamic_cast<THnSparse>(fOutput->At(i));
	207	if(hn)hn->Sumw2();
	208	}
	209
	210	TH1::AddDirectory(oldStatus);
	211
	212	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
	213	}
	214
	215	//________________________________________________________________________
	216	Bool_t AliAnalysisTaskEmcalJetTagger::Run()
	217	{
	218	// Run analysis code here, if needed. It will be executed before FillHistograms().
396ffc82	219
ba7663ad	220	MatchJetsGeo(fContainerBase,fContainerTag,0,0.3,2);
	221
	222	// if(fJetTaggingMethod==kFraction)
8612dfc8	223
	224	return kTRUE;
	225	}
	226
	227	//________________________________________________________________________
	228	Bool_t AliAnalysisTaskEmcalJetTagger::FillHistograms()
	229	{
	230	// Fill histograms.
	231
ba7663ad	232	AliEmcalJet *jet1 = NULL;
	233	AliJetContainer *jetCont = GetJetContainer(fContainerBase);
	234	if(!jetCont) return kFALSE;
	235	jetCont->ResetCurrentID();
	236	while((jet1 = jetCont->GetNextAcceptJet())) {
	237	Double_t ptJet1 = jet1->Pt() - jetCont->GetRhoVal()*jet1->Area();
8612dfc8	238	fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
8612dfc8	239
396ffc82 MV	240	//fill histo with angle between jet axis and constituents
396ffc82 MV	241	for(Int_t icc=0; icc<jet1->GetNumberOfTracks(); icc++) {
ba7663ad	242	AliVParticle vp = static_cast<AliVParticle>(jet1->TrackAt(icc, jetCont->GetParticleContainer()->GetArray()));//fTracks));
396ffc82 MV	243	if(!vp) continue;
	244	Double_t dEta = jet1->Eta()-vp->Eta();
	245	Double_t dPhi = jet1->Phi()-vp->Phi();
	246	if(dPhi<TMath::Pi()) dPhi+=TMath::TwoPi();
	247	if(dPhi>TMath::Pi()) dPhi-=TMath::TwoPi();
	248	fh3PtJetDEtaDPhiConst->Fill(ptJet1,dEta,dPhi);
	249
	250	Double_t dR = TMath::Sqrt(dPhidPhi+dEtadEta);
	251	fh2PtJetDRConst->Fill(ptJet1,dR);
	252	fh3PtJetAreaDRConst->Fill(ptJet1,jet1->Area(),dR);
	253	}
	254
8612dfc8	255	if(jet1->GetTagStatus()<1 && fJetTaggingType==kTag)
	256	continue;
	257
ba7663ad	258	AliEmcalJet *jet2 = NULL;
	259	if(fJetTaggingType==kTag) jet2 = jet1->GetTaggedJet();
	260	if(fJetTaggingType==kClosest) jet2 = jet1->ClosestJet();
8612dfc8	261	if(!jet2) continue;
ba7663ad	262
8612dfc8	263	Double_t ptJet2 = jet2->Pt() - GetRhoVal(fContainerTag)*jet2->Area();
ba7663ad	264	Double_t fraction = jetCont->GetFractionSharedPt(jet1);
ba7663ad	265	fh2PtJet2VsFraction[fCentBin]->Fill(ptJet2,fraction);
8612dfc8	266
ba7663ad	267	if(fraction<fMinFractionShared && fJetTaggingType==kClosest)
	268	continue;
	269
	270	fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
8612dfc8	271	fh2PtJet1VsPtJet2[fCentBin]->Fill(ptJet1,ptJet2);
	272
	273	Double_t dPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
	274	if(dPhi>TMath::Pi())
	275	dPhi -= TMath::TwoPi();
	276	if(dPhi<(-1.*TMath::Pi()))
	277	dPhi += TMath::TwoPi();
	278
	279	fh3PtJet1VsDeltaEtaDeltaPhi[fCentBin]->Fill(ptJet1,jet1->Eta()-jet2->Eta(),dPhi);
	280	fh2PtJet1VsDeltaR[fCentBin]->Fill(ptJet1,GetDeltaR(jet1,jet2));
	281	}
8612dfc8	282	return kTRUE;
8612dfc8	283	}
	284
	285	//________________________________________________________________________
	286	void AliAnalysisTaskEmcalJetTagger::ResetTagging(const Int_t c) {
	287
	288	//Reset tagging of container c
	289
	290	for(int i = 0;i<GetNJets(c);i++){
	291	AliEmcalJet jet = static_cast<AliEmcalJet>(GetJetFromArray(i, c));
	292	if(fJetTaggingType==kClosest)
	293	jet->ResetMatching();
	294	else if(fJetTaggingType==kTag) {
	295	jet->SetTaggedJet(0x0);
	296	jet->SetTagStatus(-1);
	297	}
	298	}
	299	}
	300
	301	//________________________________________________________________________
	302	void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
ba7663ad	303	Int_t iDebug, Float_t maxDist, Int_t type, Bool_t bReset) {
8612dfc8	304
	305	//
	306	// Match the full jets to the corresponding charged jets
	307	// Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
	308	// type: 0 = use acceptance cuts of container 1 = allow 0.1 one more for c2 in eta 2 = allow 0.1 more in eta and phi for c2
	309
	310	if(c1<0) c1 = fContainerBase;
	311	if(c2<0) c2 = fContainerTag;
	312
	313	const Int_t nJets1 = GetNJets(c1);
	314	const Int_t nJets2 = GetNJets(c2);
	315
	316	if(nJets1==0 \|\| nJets2==0) return;
	317
ba7663ad	318	if(bReset) {
	319	ResetTagging(c1);
	320	ResetTagging(c2);
	321	}
8612dfc8	322	fMatchingDone = kFALSE;
	323
	324	TArrayI faMatchIndex1;
	325	faMatchIndex1.Set(nJets2+1);
	326	faMatchIndex1.Reset(-1);
	327
	328	TArrayI faMatchIndex2;
	329	faMatchIndex2.Set(nJets1+1);
	330	faMatchIndex2.Reset(-1);
	331
	332	static TArrayS iFlag(nJets1*nJets2);
	333	if(iFlag.GetSize()<(nJets1*nJets2)){
66a08338	334	iFlag.Set(nJets1*nJets2+1);
8612dfc8	335	}
	336	iFlag.Reset(0);
	337
66a08338	338
8612dfc8	339	AliJetContainer *cont2 = GetJetContainer(c2);
	340
	341	// find the closest distance to the full jet
	342	for(int i = 0;i<nJets1;i++){
	343
	344	AliEmcalJet jet1 = static_cast<AliEmcalJet>(GetAcceptJetFromArray(i, c1));
	345	if(!jet1) continue;
	346
	347	Float_t dist = maxDist;
	348
	349	for(int j = 0;j <nJets2; j++){
	350	AliEmcalJet *jet2 = 0x0;
	351	if(type==0)
	352	jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
	353	else {
	354	jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
	355	if(!jet2) continue;
	356	if(type>0) {
	357	if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) \|\| jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
	358	continue;
	359	if(type==2) {
	360	if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) \|\| jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
	361	continue;
	362	}
	363	}
	364	}
	365	if(!jet2)
	366	continue;
	367
	368	Double_t dR = GetDeltaR(jet1,jet2);
	369	if(dR<dist && dR<maxDist){
	370	faMatchIndex2[i]=j;
	371	dist = dR;
	372	}
	373	}//j jet loop
f2baf8e8	374	if(faMatchIndex2[i]>=0) {
	375	iFlag[i*nJets2+faMatchIndex2[i]]+=1;//j closest to i
	376	if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",i,faMatchIndex2[i],dist,inJets2+faMatchIndex2[i],iFlag[inJets2+faMatchIndex2[i]]);
	377	}
8612dfc8	378
	379	}//i jet loop
	380
	381
	382	// other way around
	383	for(int j = 0;j<nJets2;j++){
	384	AliEmcalJet *jet2 = 0x0;
	385	if(type==0)
	386	jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
	387	else {
	388	jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
	389	if(!jet2) continue;;
	390	if(type>0) {
	391	if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) \|\| jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
	392	continue;
	393	if(type==2) {
	394	if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) \|\| jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
	395	continue;
	396	}
	397	}
	398	}
	399	if(!jet2)
	400	continue;
	401
	402	Float_t dist = maxDist;
	403	for(int i = 0;i<nJets1;i++){
	404	AliEmcalJet jet1 = static_cast<AliEmcalJet>(GetAcceptJetFromArray(i, c1));
	405	if(!jet1)
	406	continue;
	407
	408	Double_t dR = GetDeltaR(jet1,jet2);
	409	if(dR<dist && dR<maxDist){
	410	faMatchIndex1[j]=i;
	411	dist = dR;
	412	}
	413	}
f2baf8e8	414	if(faMatchIndex1[j]>=0) {
	415	iFlag[faMatchIndex1[j]*nJets2+j]+=2;//i closest to j
	416	if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faMatchIndex1[j],j,dist,faMatchIndex1[j]nJets2+j,iFlag[faMatchIndex1[j]nJets2+j]);
	417	}
8612dfc8	418
	419	}
	420
	421	// check for "true" correlations
	422	for(int i = 0;i<nJets1;i++){
	423	AliEmcalJet jet1 = static_cast<AliEmcalJet>(GetJetFromArray(i, c1));
	424	for(int j = 0;j<nJets2;j++){
	425	AliEmcalJet jet2 = static_cast<AliEmcalJet>(GetJetFromArray(j, c2));
f2baf8e8	426	AliDebug(11,Form("%s: Flag[%d][%d] %d ",GetName(),i,j,iFlag[i*nJets2+j]));
8612dfc8	427
	428	// we have a uniqe correlation
	429	if(iFlag[i*nJets2+j]==3){
	430
	431	Double_t dR = GetDeltaR(jet1,jet2);
	432	if(iDebug>1) Printf("closest jets %d %d dR = %f",j,i,dR);
	433
	434	if(fJetTaggingType==kTag) {
	435	jet1->SetTaggedJet(jet2);
	436	jet1->SetTagStatus(1);
	437
	438	jet2->SetTaggedJet(jet1);
	439	jet2->SetTagStatus(1);
	440	}
	441	else if(fJetTaggingType==kClosest) {
	442	jet1->SetClosestJet(jet2,dR);
	443	jet2->SetClosestJet(jet1,dR);
	444	}
	445	}
	446	}
	447	}
8612dfc8	448	fMatchingDone = kTRUE;
8612dfc8	449	}
	450
	451	//________________________________________________________________________
	452	Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
	453	//
	454	// Helper function to calculate the distance between two jets
	455	//
	456
	457	Double_t dPhi = jet1->Phi() - jet2->Phi();
	458	if(dPhi>TMath::Pi())
	459	dPhi -= TMath::TwoPi();
	460	if(dPhi<(-1.*TMath::Pi()))
	461	dPhi += TMath::TwoPi();
	462	Double_t dEta = jet1->Eta() - jet2->Eta();
	463	Double_t dR = TMath::Sqrt(dPhidPhi+dEtadEta);
	464	return dR;
	465	}
	466
	467	//________________________________________________________________________
	468	Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
	469	//
	470	// Calculate azimuthal angle between the axises of the jets
	471	//
	472
	473	return GetDeltaPhi(jet1->Phi(),jet2->Phi());
	474
	475	}
	476
	477	//________________________________________________________________________
	478	Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(Double_t phi1,Double_t phi2) {
	479	//
	480	// Calculate azimuthal angle between the axises of the jets
	481	//
	482
	483	Double_t dPhi = phi1-phi2;
	484
	485	if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
	486	if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
	487
	488	return dPhi;
	489	}
	490
	491	//________________________________________________________________________
	492	Double_t AliAnalysisTaskEmcalJetTagger::GetFractionSharedPt(const AliEmcalJet jet1, const AliEmcalJet jet2) const
	493	{
	494	//
	495	// Get fraction of shared pT between matched full and charged jet
	496	// Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
	497	//
	498
	499	Double_t fraction = 0.;
	500	Double_t jetPt2 = jet2->Pt();
	501
	502	if(jetPt2>0) {
	503
ba7663ad	504	AliDebug(2,Form("%s: nConstituents_jet1: %d, nConstituents_jet2: %d, tracks_jet1: %d tracks_jet2: %d clusters_jet1: %d clusters_jet2: %d",GetName(),jet1->GetNumberOfConstituents(),jet2->GetNumberOfConstituents(),jet1->GetNumberOfTracks(),jet2->GetNumberOfTracks(),jet1->GetNumberOfClusters(),jet2->GetNumberOfClusters()));
8612dfc8	505
	506	Double_t sumPt = 0.;
	507	AliVParticle *vpf = 0x0;
	508	Int_t iFound = 0;
	509	for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
	510	Int_t idx = (Int_t)jet2->TrackAt(icc);
	511	iFound = 0;
	512	for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
	513	if(idx == jet1->TrackAt(icf) && iFound==0 ) {
	514	iFound=1;
	515	vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fTracks));
ba7663ad	516	if(vpf) sumPt += vpf->Pt();
8612dfc8	517	continue;
	518	}
	519	}
	520	}
8612dfc8	521	fraction = sumPt/jetPt2;
ba7663ad	522	} else fraction = -1.;
8612dfc8	523	return fraction;
8612dfc8	524	}
	525
	526	//________________________________________________________________________
	527	Bool_t AliAnalysisTaskEmcalJetTagger::RetrieveEventObjects() {
	528	//
	529	// retrieve event objects
	530	//
	531
	532	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
	533	return kFALSE;
	534
	535	return kTRUE;
	536
	537	}
	538
	539	//_______________________________________________________________________
	540	void AliAnalysisTaskEmcalJetTagger::Terminate(Option_t *)
	541	{
	542	// Called once at the end of the analysis.
	543	}
	544