[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),
  fh3PtJet1VsMassJet1MassJet2(0),
  fh3PtJet2VsMassJet1MassJet2(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];
  fh3PtJet1VsMassJet1MassJet2  = new TH3F*[fNcentBins];
  fh3PtJet2VsMassJet1MassJet2  = new TH3F*[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;
    fh3PtJet1VsMassJet1MassJet2[i] = 0;
    fh3PtJet2VsMassJet1MassJet2[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),
  fh3PtJet1VsMassJet1MassJet2(0),
  fh3PtJet2VsMassJet1MassJet2(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];
  fh3PtJet1VsMassJet1MassJet2 = new TH3F*[fNcentBins];
  fh3PtJet2VsMassJet1MassJet2 = new TH3F*[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;
    fh3PtJet1VsMassJet1MassJet2[i] = 0;
    fh3PtJet2VsMassJet1MassJet2[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 Int_t nBinsM           = 150;

  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;
  const Double_t minM        = -20.;
  const Double_t maxM        = 130.;

  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]);

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

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