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

//
// Dijet base 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 "AliAnalysisTaskEmcalDiJetBase.h"

ClassImp(AliAnalysisTaskEmcalDiJetBase)

//________________________________________________________________________
AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalDiJetBase", kTRUE),
  fDebug(kFALSE),
  fJetCorrelationType(kCorrelateAll),
  fJetFullChargedMatchingType(kFraction),
  fTriggerClass(""),
  fContainerCharged(1),
  fContainerFull(0),
  fContainerChargedMC(3),
  fContainerFullMC(2),
  fRhoType(0),
  fRhoChVal(0),
  fRhoFullVal(0),
  fDoChargedCharged(kTRUE),
  fDoFullCharged(kTRUE),
  fDoFullFull(kFALSE),
  fPtMinTriggerJet(10.),
  fMinFractionShared(0.5),
  fMatchingDone(kFALSE),
  faFullFracIndex(0),
  faFullFracIndexMC(0),
  fhNEvents(0),
  fHistTrialsSelEvents(0)
{
  // Default constructor.

  
  SetMakeGeneralHistograms(kTRUE);

}

//________________________________________________________________________
AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),
  fDebug(kFALSE),
  fJetCorrelationType(kCorrelateAll),
  fJetFullChargedMatchingType(kFraction),
  fTriggerClass(""),
  fContainerCharged(1),
  fContainerFull(0),
  fContainerChargedMC(3),
  fContainerFullMC(2),
  fRhoType(0),
  fRhoChVal(0),
  fRhoFullVal(0),
  fDoChargedCharged(kTRUE),
  fDoFullCharged(kTRUE),
  fDoFullFull(kFALSE),
  fPtMinTriggerJet(10.),
  fMinFractionShared(0.5),
  fMatchingDone(kFALSE),
  faFullFracIndex(0),
  faFullFracIndexMC(0),
  fhNEvents(0),
  fHistTrialsSelEvents(0)
{
  // Standard constructor.

  SetMakeGeneralHistograms(kTRUE);

}

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

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalDiJetBase::SelectEvent() {
  //
  // Decide if event should be selected for analysis
  //

  fhNEvents->Fill(3.5);

  if(!fTriggerClass.IsNull()) {
    //Check if requested trigger was fired
    TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
    if(!firedTrigClass.Contains(fTriggerClass))
      return kFALSE;
    if(fTriggerClass.Contains("J1") && firedTrigClass.Contains("J2")) //assign J1&J2 triggers to J2 category
      return kFALSE;
  }

  fhNEvents->Fill(1.5);

  fHistTrialsSelEvents->Fill(fPtHardBin, fNTrials);
  
  return kTRUE;

}

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

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

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

  fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
  fOutput->Add(fhNEvents);

  //Pythia info
  fHistTrialsSelEvents = new TH1F("fHistTrialsSelEvents", "fHistTrialsSelEvents", 12, 0, 12);
  fHistTrialsSelEvents->GetXaxis()->SetTitle("p_{T} hard bin");
  fHistTrialsSelEvents->GetYaxis()->SetTitle("trials");
  fOutput->Add(fHistTrialsSelEvents);

  TH1::AddDirectory(oldStatus);

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

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalDiJetBase::IsSameJet(const Int_t ijt, const Int_t ija, const Int_t type, const Bool_t isMC) {
   //check if two jets are the same one

   Bool_t bSame = kFALSE;

   if(type<2 && ijt==ija)
     bSame = kTRUE;
   if(type==2) {

     if(fJetFullChargedMatchingType==kFraction) {
       if(isMC && faFullFracIndexMC.At(ijt)==ija)
	 bSame = kTRUE;
       else if(!isMC && faFullFracIndex.At(ijt)==ija)
	 bSame = kTRUE;
     }
     else if(fJetFullChargedMatchingType==kGeo) {
       Int_t contFull = fContainerFull;
       Int_t contChar = fContainerCharged;
       if(isMC) {
	 contFull = fContainerFullMC;
	 contChar = fContainerChargedMC;
       }
       AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, contFull));
       AliEmcalJet *chJet   = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, contChar));
       AliEmcalJet *matchJet = fullJet->ClosestJet();
       if(chJet==matchJet)
	 bSame = kTRUE;
     }
     else if(fJetFullChargedMatchingType==kNoMatching) {
       return kFALSE;
     }
     else{
       AliWarning(Form("%s: matching type unknown", GetName()));
     }
   }

   return bSame;
 }


//________________________________________________________________________
Double_t AliAnalysisTaskEmcalDiJetBase::GetJetPt(const AliEmcalJet *jet, const Int_t type) {

  if(!jet) return -99;

  if(type==0)
    return jet->Pt() - fRhoFullVal * jet->Area();
  else if(type==1)
    return jet->Pt() - fRhoChVal * jet->Area();
  else
    return jet->Pt();
  
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const AliVParticle *trk, const AliEmcalJet *jet)          const
{  
  // Get Z of constituent trk

  return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
{
  // 
  // Get the z of a constituent inside of a jet
  //

  Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;

  if(pJetSq>0.)
    return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
  else {
    AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq)); 
    return 0;
  }
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalDiJetBase::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 = dPhi - 2.*TMath::Pi();
  if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
  Double_t dEta = jet1->Eta() - jet2->Eta();
  Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
  return dR;
}

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

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

}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalDiJetBase::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 AliAnalysisTaskEmcalDiJetBase::GetFractionSharedPt(const AliEmcalJet *jetFull, const AliEmcalJet *jetCharged) 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 jetPtCh = jetCharged->Pt();
 
  if(jetPtCh>0) {

    if(fDebug>10)  AliInfo(Form("%s: nConstituents: %d, ch: %d  chne: %d ne: %d",GetName(),jetFull->GetNumberOfConstituents(),jetCharged->GetNumberOfTracks(),jetFull->GetNumberOfTracks(),jetFull->GetNumberOfClusters()));
    
    Double_t sumPt = 0.;
    AliVParticle *vpf = 0x0;
    Int_t iFound = 0;
    for(Int_t icc=0; icc<jetCharged->GetNumberOfTracks(); icc++) {
      Int_t idx = (Int_t)jetCharged->TrackAt(icc);
      iFound = 0;
      for(Int_t icf=0; icf<jetFull->GetNumberOfTracks(); icf++) {
	if(idx == jetFull->TrackAt(icf) && iFound==0 ) {
	  iFound=1;
	  vpf = static_cast<AliVParticle*>(jetFull->TrackAt(icf, fTracks));
	  if(!vpf) continue;
	  if(vpf->Charge()!=0)
	    sumPt += vpf->Pt();
	  continue;
	}
      }
    }
    
    fraction = sumPt/jetPtCh;
  }

  if(fDebug>10) AliInfo(Form("%s: charged shared fraction: %.2f",GetName(),fraction));

  return fraction;

}

//________________________________________________________________________
void AliAnalysisTaskEmcalDiJetBase::MatchJetsGeo(Int_t cFull, Int_t cCharged,
						 Int_t iDebug, Float_t maxDist, Int_t type) {

  //
  // 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 cCharged(MC) in eta 2 = allow 0.1 more in eta and phi for cCharged(MC)

  const int kMode = 3;

  const Int_t nChJets   = GetNJets(cCharged);
  const Int_t nFullJets = GetNJets(cFull);

  if(nFullJets==0 || nChJets==0) return;

  TArrayI faChNeMatchIndex;
  faChNeMatchIndex.Set(nChJets+1);
  faChNeMatchIndex.Reset(-1);

  TArrayI faChMatchIndex;
  faChMatchIndex.Set(nFullJets+1);
  faChMatchIndex.Reset(-1);

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

  AliJetContainer *contCh = GetJetContainer(cCharged);

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

    AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
    if(!fullJet) continue;

    Float_t dist = maxDist;
    
    for(int ich = 0;ich<nChJets;ich++){
      AliEmcalJet *chJet = 0x0;
      if(type==0)
	chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
      else {
	chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
	if(!chJet) continue;
	if(type>0) {
	  if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
	    continue;
	  if(type==2) {
	    if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
	      continue;
	  }
	}
      }
      if(!chJet)
	continue;

      Double_t frac = GetFractionSharedPt(fullJet,chJet);
      Double_t dR = GetDeltaR(fullJet,chJet);
      if(dR<dist){
	faChMatchIndex[ifu]=ich;
	dist = dR;
      }
      if(iDebug>10) Printf("Distance (%d)--(%d) %3.3f frac:%.2f",ifu,ich,dR,frac);
    }
    if(faChMatchIndex[ifu]>=0) iFlag[ifu*nChJets+faChMatchIndex[ifu]]+=1;//ich closest to ifu
    if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",ifu,faChMatchIndex[ifu],dist,ifu*nChJets+faChMatchIndex[ifu],iFlag[ifu*nChJets+faChMatchIndex[ifu]]);
    
    // reset...
    faChMatchIndex[ifu]=-1;

           
  }//full jet loop


  // other way around
  for(int ich = 0;ich<nChJets;ich++){
    AliEmcalJet *chJet = 0x0;
    if(type==0)
      chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
    else {
      chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
      if(!chJet) continue;;
      if(type>0) {
	if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
	  continue;
	if(type==2) {
	  if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
	    continue;
	}
      }
    }
    if(!chJet)
      continue;

    Float_t dist = maxDist;
    for(int ifu = 0;ifu<nFullJets;ifu++){
      AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
      if(!fullJet)
	continue;

      Double_t dR = GetDeltaR(fullJet,chJet); 
      if(dR<dist){
	faChNeMatchIndex[ich]=ifu;
        dist = dR;
      }   
    }
    if(faChNeMatchIndex[ich]>=0) iFlag[faChNeMatchIndex[ich]*nChJets+ich]+=2;//ifu closest to ich
    if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faChNeMatchIndex[ich],ich,dist,faChNeMatchIndex[ich]*nChJets+ich,iFlag[faChNeMatchIndex[ich]*nChJets+ich]);

    // reset...
    faChNeMatchIndex[ich]=-1;
    
  }
  
  
  // check for "true" correlations
  for(int ifu = 0;ifu<nFullJets;ifu++){
    for(int ich = 0;ich<nChJets;ich++){
      if(iDebug>10) AliInfo(Form("%s: Flag[%d][%d] %d ",GetName(),ifu,ich,iFlag[ifu*nChJets+ich]));
      
      if(kMode==3){
        // we have a uniqe correlation
        if(iFlag[ifu*nChJets+ich]==3){

	  AliEmcalJet *chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
	  AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetJetFromArray(ifu, cFull));
	  Double_t dR = GetDeltaR(fullJet,chJet); 

	  if(iDebug>10) Printf("closest jets %d  %d  dR =  %f",ich,ifu,dR);

	  chJet->SetClosestJet(fullJet,dR);
	  fullJet->SetClosestJet(chJet,dR);

	}
      }
    }
  }

  fMatchingDone = kTRUE;
  
}

//________________________________________________________________________
void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndex() {

  // take each full jet and set the index of the charged jet with the largest shared charged fraction

  const Int_t nJetsCh   = GetNJets(fContainerCharged);
  const Int_t nJetsFull = GetNJets(fContainerFull);
  faFullFracIndex.Set(nJetsFull+1);
  faFullFracIndex.Reset(-1);

  AliJetContainer  *cont = GetJetContainer(fContainerFull); 
  Float_t radius =  cont->GetJetRadius();

  for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
    Double_t frac = 0.;
    Double_t dist = 10.;
    AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFull));
    if(!jetFull) {
      faFullFracIndex.SetAt(-1,ifu);
      continue;
    }

    for(Int_t ich = 0; ich<nJetsCh; ich++) {
      AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerCharged));
      if(!jetCh)
	continue;
      Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
      dist = GetDeltaR(jetFull,jetCh);
      if(tmpFrac>frac && dist<radius) {
	frac = tmpFrac;
	faFullFracIndex.SetAt(ich,ifu);
      }
    }

  }

}

//________________________________________________________________________
void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndexMC() {

  // take each full jet and set the index of the charged jet with the largest shared charged fraction

  const Int_t nJetsCh   = GetNJets(fContainerChargedMC);
  const Int_t nJetsFull = GetNJets(fContainerFullMC);
  faFullFracIndexMC.Set(nJetsFull);
  faFullFracIndexMC.Reset(-1);

  AliJetContainer  *cont = GetJetContainer(fContainerFullMC); 
  Float_t radius =  cont->GetJetRadius();

  for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
    Double_t frac = 0.;
    Double_t dist = 10.;
    AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFullMC));
    if(!jetFull) {
      faFullFracIndexMC.SetAt(-1,ifu);
      continue;
    }

    for(Int_t ich = 0; ich<nJetsCh; ich++) {
      AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerChargedMC));
      if(!jetCh)
	continue;
      Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
      dist = GetDeltaR(jetFull,jetCh);
      if(tmpFrac>frac && dist<radius) {
	frac = tmpFrac;
	faFullFracIndexMC.SetAt(ich,ifu);
      }
    }

  }

}

//_______________________________________________________________________
AliEmcalJet* AliAnalysisTaskEmcalDiJetBase::GetLeadingJetOppositeHemisphere(const Int_t type, const Int_t typea, const AliEmcalJet *jetTrig) {

  // Get leading jet in opposite hemisphere from trigger jet
  // type = correlation type
  // typea = container of associated jets

  Int_t nJetsAssoc = GetNJets(typea);
  Double_t ptLead = -999;
  Int_t    iJetLead = -1;
  for(Int_t ija=0; ija<nJetsAssoc; ija++) {

    AliEmcalJet *jetAssoc = NULL;
    if(type==0) {
      jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
      if(TMath::Abs(jetAssoc->Eta())>0.5)
        jetAssoc = NULL;
    }
    else
      jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));

    if(!jetAssoc)
      continue;

    Double_t dPhi = GetDeltaPhi(jetTrig,jetAssoc);
    Double_t phiMin = 0.5*TMath::Pi();
    Double_t phiMax = 1.5*TMath::Pi();
    if(dPhi<phiMin || dPhi>phiMax)
      continue;

    Double_t jetAssocPt = GetJetPt(jetAssoc,typea);

    if(jetAssocPt>ptLead) {
      ptLead = jetAssocPt;
      iJetLead = ija;
    }

  }

  AliEmcalJet *jetAssocLead = NULL;
  if(iJetLead>-1)
    jetAssocLead = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead, typea));

  return jetAssocLead;

}

//_______________________________________________________________________
AliEmcalJet* AliAnalysisTaskEmcalDiJetBase::GetSecondLeadingJetOppositeHemisphere(const Int_t type, const Int_t typea, const AliEmcalJet *jetTrig) {

  // Get leading jet in opposite hemisphere from trigger jet
  // type = correlation type
  // typea = container of associated jets

  Int_t nJetsAssoc = GetNJets(typea);
  Double_t ptLead = -999;
  Int_t    iJetLead = -1;
  Int_t    iJetLead2 = -1;
  for(Int_t ija=0; ija<nJetsAssoc; ija++) {

    AliEmcalJet *jetAssoc = NULL;
    if(type==0) {
      jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
      if(TMath::Abs(jetAssoc->Eta())>0.5)
        jetAssoc = NULL;
    }
    else
      jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));

    if(!jetAssoc)
      continue;

    Double_t dPhi = GetDeltaPhi(jetTrig,jetAssoc);
    Double_t phiMin = 0.5*TMath::Pi();
    Double_t phiMax = 1.5*TMath::Pi();
    if(dPhi<phiMin || dPhi>phiMax)
      continue;

    Double_t jetAssocPt = GetJetPt(jetAssoc,typea);

    if(jetAssocPt>ptLead) {
      iJetLead2 = iJetLead;
      ptLead = jetAssocPt;
      iJetLead = ija;
    }

  }

  AliEmcalJet *jetAssocLead2 = NULL;
  if(iJetLead2>-1)
    jetAssocLead2 = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead2, typea));

  return jetAssocLead2;

}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects() {
  //
  // Retrieve objects from event.
  //

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

  if(fRhoType==0) {
    fRhoFullVal = 0.;
    fRhoChVal = 0.;
  }
  if(fRhoType==1) {
    fRhoFullVal = GetRhoVal(fContainerFull);
    fRhoChVal = GetRhoVal(fContainerCharged);
  }

  return kTRUE;
}

//_______________________________________________________________________
void AliAnalysisTaskEmcalDiJetBase::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
6e8b6371	1	//
	2	// Dijet base 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"
6e8b6371	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"
	31
	32	#include "AliAnalysisTaskEmcalDiJetBase.h"
	33
	34	ClassImp(AliAnalysisTaskEmcalDiJetBase)
	35
	36	//________________________________________________________________________
	37	AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase() :
9239b066	38	AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalDiJetBase", kTRUE),
6e8b6371	39	fDebug(kFALSE),
c66364a6	40	fJetCorrelationType(kCorrelateAll),
6e8b6371	41	fJetFullChargedMatchingType(kFraction),
	42	fTriggerClass(""),
	43	fContainerCharged(1),
	44	fContainerFull(0),
	45	fContainerChargedMC(3),
	46	fContainerFullMC(2),
	47	fRhoType(0),
	48	fRhoChVal(0),
	49	fRhoFullVal(0),
	50	fDoChargedCharged(kTRUE),
	51	fDoFullCharged(kTRUE),
bb84b374	52	fDoFullFull(kFALSE),
6e8b6371	53	fPtMinTriggerJet(10.),
	54	fMinFractionShared(0.5),
	55	fMatchingDone(kFALSE),
	56	faFullFracIndex(0),
	57	faFullFracIndexMC(0),
6e8b6371	58	fhNEvents(0),
8e49a788	59	fHistTrialsSelEvents(0)
6e8b6371	60	{
	61	// Default constructor.
	62
	63
	64	SetMakeGeneralHistograms(kTRUE);
	65
	66	}
	67
	68	//________________________________________________________________________
	69	AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase(const char *name) :
9239b066	70	AliAnalysisTaskEmcalJet(name, kTRUE),
6e8b6371	71	fDebug(kFALSE),
c66364a6	72	fJetCorrelationType(kCorrelateAll),
6e8b6371	73	fJetFullChargedMatchingType(kFraction),
	74	fTriggerClass(""),
	75	fContainerCharged(1),
	76	fContainerFull(0),
	77	fContainerChargedMC(3),
	78	fContainerFullMC(2),
	79	fRhoType(0),
	80	fRhoChVal(0),
	81	fRhoFullVal(0),
	82	fDoChargedCharged(kTRUE),
	83	fDoFullCharged(kTRUE),
bb84b374	84	fDoFullFull(kFALSE),
6e8b6371	85	fPtMinTriggerJet(10.),
	86	fMinFractionShared(0.5),
	87	fMatchingDone(kFALSE),
	88	faFullFracIndex(0),
	89	faFullFracIndexMC(0),
6e8b6371	90	fhNEvents(0),
8e49a788	91	fHistTrialsSelEvents(0)
6e8b6371	92	{
	93	// Standard constructor.
	94
	95	SetMakeGeneralHistograms(kTRUE);
	96
	97	}
	98
	99	//________________________________________________________________________
	100	AliAnalysisTaskEmcalDiJetBase::~AliAnalysisTaskEmcalDiJetBase()
	101	{
	102	// Destructor.
	103	}
	104
6e8b6371	105	//________________________________________________________________________
	106	Bool_t AliAnalysisTaskEmcalDiJetBase::SelectEvent() {
	107	//
	108	// Decide if event should be selected for analysis
	109	//
	110
6e8b6371	111	fhNEvents->Fill(3.5);
	112
	113	if(!fTriggerClass.IsNull()) {
	114	//Check if requested trigger was fired
	115	TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
	116	if(!firedTrigClass.Contains(fTriggerClass))
	117	return kFALSE;
a29bc208	118	if(fTriggerClass.Contains("J1") && firedTrigClass.Contains("J2")) //assign J1&J2 triggers to J2 category
8e49a788	119	return kFALSE;
6e8b6371	120	}
6e8b6371	121
a29bc208	122	fhNEvents->Fill(1.5);
8e49a788	123
8e49a788	124	fHistTrialsSelEvents->Fill(fPtHardBin, fNTrials);
6e8b6371	125
	126	return kTRUE;
	127
	128	}
	129
	130	//________________________________________________________________________
	131	void AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects()
	132	{
	133	// Create user output.
	134
9239b066	135	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
6e8b6371	136
	137	Bool_t oldStatus = TH1::AddDirectoryStatus();
	138	TH1::AddDirectory(kFALSE);
	139
	140	fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
	141	fOutput->Add(fhNEvents);
	142
6e8b6371	143	//Pythia info
6e8b6371	144	fHistTrialsSelEvents = new TH1F("fHistTrialsSelEvents", "fHistTrialsSelEvents", 12, 0, 12);
	145	fHistTrialsSelEvents->GetXaxis()->SetTitle("p_{T} hard bin");
	146	fHistTrialsSelEvents->GetYaxis()->SetTitle("trials");
	147	fOutput->Add(fHistTrialsSelEvents);
6e8b6371	148
	149	TH1::AddDirectory(oldStatus);
	150
	151	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
	152	}
	153
	154	//________________________________________________________________________
	155	Bool_t AliAnalysisTaskEmcalDiJetBase::IsSameJet(const Int_t ijt, const Int_t ija, const Int_t type, const Bool_t isMC) {
	156	//check if two jets are the same one
	157
	158	Bool_t bSame = kFALSE;
	159
	160	if(type<2 && ijt==ija)
	161	bSame = kTRUE;
	162	if(type==2) {
	163
	164	if(fJetFullChargedMatchingType==kFraction) {
	165	if(isMC && faFullFracIndexMC.At(ijt)==ija)
	166	bSame = kTRUE;
	167	else if(!isMC && faFullFracIndex.At(ijt)==ija)
	168	bSame = kTRUE;
	169	}
	170	else if(fJetFullChargedMatchingType==kGeo) {
	171	Int_t contFull = fContainerFull;
	172	Int_t contChar = fContainerCharged;
	173	if(isMC) {
	174	contFull = fContainerFullMC;
	175	contChar = fContainerChargedMC;
	176	}
	177	AliEmcalJet fullJet = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ijt, contFull));
	178	AliEmcalJet chJet = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ija, contChar));
	179	AliEmcalJet *matchJet = fullJet->ClosestJet();
	180	if(chJet==matchJet)
	181	bSame = kTRUE;
	182	}
	183	else if(fJetFullChargedMatchingType==kNoMatching) {
	184	return kFALSE;
	185	}
	186	else{
	187	AliWarning(Form("%s: matching type unknown", GetName()));
	188	}
6e8b6371	189	}
	190
	191	return bSame;
	192	}
	193
	194
	195	//________________________________________________________________________
	196	Double_t AliAnalysisTaskEmcalDiJetBase::GetJetPt(const AliEmcalJet *jet, const Int_t type) {
	197
	198	if(!jet) return -99;
	199
	200	if(type==0)
	201	return jet->Pt() - fRhoFullVal * jet->Area();
	202	else if(type==1)
	203	return jet->Pt() - fRhoChVal * jet->Area();
	204	else
	205	return jet->Pt();
	206
	207	}
	208
	209	//________________________________________________________________________
	210	Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const AliVParticle trk, const AliEmcalJet jet) const
	211	{
	212	// Get Z of constituent trk
	213
	214	return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
	215	}
	216
	217	//________________________________________________________________________
	218	Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
	219	{
	220	//
	221	// Get the z of a constituent inside of a jet
	222	//
	223
	224	Double_t pJetSq = jetPxjetPx+jetPyjetPy+jetPz*jetPz;
	225
	226	if(pJetSq>0.)
	227	return (trkPxjetPx+trkPyjetPy+trkPz*jetPz)/pJetSq;
	228	else {
	229	AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
	230	return 0;
	231	}
	232	}
	233
	234	//________________________________________________________________________
	235	Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
	236	//
	237	// Helper function to calculate the distance between two jets
	238	//
	239
	240	Double_t dPhi = jet1->Phi() - jet2->Phi();
	241	if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
	242	if(dPhi<(-1.TMath::Pi()))dPhi = dPhi + 2.TMath::Pi();
	243	Double_t dEta = jet1->Eta() - jet2->Eta();
	244	Double_t dR = TMath::Sqrt(dPhidPhi+dEtadEta);
	245	return dR;
	246	}
	247
	248	//________________________________________________________________________
	249	Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
	250	//
	251	// Calculate azimuthal angle between the axises of the jets
	252	//
253
254	return GetDeltaPhi(jet1->Phi(),jet2->Phi());
255
256	}
257
258	//________________________________________________________________________
259	Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(Double_t phi1,Double_t phi2) {
260	//
261	// Calculate azimuthal angle between the axises of the jets
262	//
263
264	Double_t dPhi = phi1-phi2;
265
6e8b6371	266	if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
	267	if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
	268
	269	return dPhi;
	270	}
	271
	272	//________________________________________________________________________
	273	Double_t AliAnalysisTaskEmcalDiJetBase::GetFractionSharedPt(const AliEmcalJet jetFull, const AliEmcalJet jetCharged) const
	274	{
	275	//
	276	// Get fraction of shared pT between matched full and charged jet
	277	// Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
	278	//
	279
	280	Double_t fraction = 0.;
	281	Double_t jetPtCh = jetCharged->Pt();
	282
	283	if(jetPtCh>0) {
	284
	285	if(fDebug>10) AliInfo(Form("%s: nConstituents: %d, ch: %d chne: %d ne: %d",GetName(),jetFull->GetNumberOfConstituents(),jetCharged->GetNumberOfTracks(),jetFull->GetNumberOfTracks(),jetFull->GetNumberOfClusters()));
	286
	287	Double_t sumPt = 0.;
	288	AliVParticle *vpf = 0x0;
	289	Int_t iFound = 0;
	290	for(Int_t icc=0; icc<jetCharged->GetNumberOfTracks(); icc++) {
	291	Int_t idx = (Int_t)jetCharged->TrackAt(icc);
	292	iFound = 0;
	293	for(Int_t icf=0; icf<jetFull->GetNumberOfTracks(); icf++) {
	294	if(idx == jetFull->TrackAt(icf) && iFound==0 ) {
	295	iFound=1;
	296	vpf = static_cast<AliVParticle*>(jetFull->TrackAt(icf, fTracks));
	297	if(!vpf) continue;
	298	if(vpf->Charge()!=0)
	299	sumPt += vpf->Pt();
	300	continue;
	301	}
	302	}
	303	}
	304
	305	fraction = sumPt/jetPtCh;
	306	}
	307
	308	if(fDebug>10) AliInfo(Form("%s: charged shared fraction: %.2f",GetName(),fraction));
	309
	310	return fraction;
	311
	312	}
	313
	314	//________________________________________________________________________
	315	void AliAnalysisTaskEmcalDiJetBase::MatchJetsGeo(Int_t cFull, Int_t cCharged,
	316	Int_t iDebug, Float_t maxDist, Int_t type) {
	317
	318	//
	319	// Match the full jets to the corresponding charged jets
	320	// Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
	321	// type: 0 = use acceptance cuts of container 1 = allow 0.1 one more for cCharged(MC) in eta 2 = allow 0.1 more in eta and phi for cCharged(MC)
	322
	323	const int kMode = 3;
	324
	325	const Int_t nChJets = GetNJets(cCharged);
	326	const Int_t nFullJets = GetNJets(cFull);
	327
	328	if(nFullJets==0 \|\| nChJets==0) return;
	329
330	TArrayI faChNeMatchIndex;
331	faChNeMatchIndex.Set(nChJets+1);
254370e1	332	faChNeMatchIndex.Reset(-1);
6e8b6371	333
	334	TArrayI faChMatchIndex;
	335	faChMatchIndex.Set(nFullJets+1);
254370e1	336	faChMatchIndex.Reset(-1);
6e8b6371	337
	338	static TArrayS iFlag(nChJets*nFullJets);
	339	if(iFlag.GetSize()<(nChJets*nFullJets)){
	340	iFlag.Set(nChJets*nFullJets+1);
	341	}
	342	iFlag.Reset(0);
	343
	344	AliJetContainer *contCh = GetJetContainer(cCharged);
	345
	346	// find the closest distance to the full jet
	347	for(int ifu = 0;ifu<nFullJets;ifu++){
	348
	349	AliEmcalJet fullJet = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ifu, cFull));
	350	if(!fullJet) continue;
	351
	352	Float_t dist = maxDist;
	353
	354	for(int ich = 0;ich<nChJets;ich++){
	355	AliEmcalJet *chJet = 0x0;
	356	if(type==0)
	357	chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
	358	else {
	359	chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
254370e1	360	if(!chJet) continue;
6e8b6371	361	if(type>0) {
	362	if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) \|\| chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
	363	continue;
	364	if(type==2) {
	365	if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) \|\| chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
	366	continue;
	367	}
	368	}
	369	}
	370	if(!chJet)
	371	continue;
	372
	373	Double_t frac = GetFractionSharedPt(fullJet,chJet);
	374	Double_t dR = GetDeltaR(fullJet,chJet);
	375	if(dR<dist){
	376	faChMatchIndex[ifu]=ich;
	377	dist = dR;
	378	}
	379	if(iDebug>10) Printf("Distance (%d)--(%d) %3.3f frac:%.2f",ifu,ich,dR,frac);
	380	}
	381	if(faChMatchIndex[ifu]>=0) iFlag[ifu*nChJets+faChMatchIndex[ifu]]+=1;//ich closest to ifu
	382	if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",ifu,faChMatchIndex[ifu],dist,ifunChJets+faChMatchIndex[ifu],iFlag[ifunChJets+faChMatchIndex[ifu]]);
	383
	384	// reset...
	385	faChMatchIndex[ifu]=-1;
	386
	387
	388	}//full jet loop
	389
	390
	391	// other way around
	392	for(int ich = 0;ich<nChJets;ich++){
	393	AliEmcalJet *chJet = 0x0;
	394	if(type==0)
	395	chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
	396	else {
	397	chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
	398	if(!chJet) continue;;
	399	if(type>0) {
	400	if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) \|\| chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
	401	continue;
	402	if(type==2) {
	403	if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) \|\| chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
	404	continue;
	405	}
	406	}
	407	}
	408	if(!chJet)
	409	continue;
	410
	411	Float_t dist = maxDist;
	412	for(int ifu = 0;ifu<nFullJets;ifu++){
	413	AliEmcalJet fullJet = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ifu, cFull));
	414	if(!fullJet)
	415	continue;
	416
	417	Double_t dR = GetDeltaR(fullJet,chJet);
	418	if(dR<dist){
	419	faChNeMatchIndex[ich]=ifu;
	420	dist = dR;
	421	}
	422	}
	423	if(faChNeMatchIndex[ich]>=0) iFlag[faChNeMatchIndex[ich]*nChJets+ich]+=2;//ifu closest to ich
	424	if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faChNeMatchIndex[ich],ich,dist,faChNeMatchIndex[ich]nChJets+ich,iFlag[faChNeMatchIndex[ich]nChJets+ich]);
425
426	// reset...
427	faChNeMatchIndex[ich]=-1;
428
429	}
430
431
432	// check for "true" correlations
433	for(int ifu = 0;ifu<nFullJets;ifu++){
434	for(int ich = 0;ich<nChJets;ich++){
435	if(iDebug>10) AliInfo(Form("%s: Flag[%d][%d] %d ",GetName(),ifu,ich,iFlag[ifu*nChJets+ich]));
436
437	if(kMode==3){
438	// we have a uniqe correlation
439	if(iFlag[ifu*nChJets+ich]==3){
440
441	AliEmcalJet chJet = static_cast<AliEmcalJet>(GetJetFromArray(ich, cCharged));
442	AliEmcalJet fullJet = static_cast<AliEmcalJet>(GetJetFromArray(ifu, cFull));
443	Double_t dR = GetDeltaR(fullJet,chJet);
444
445	if(iDebug>10) Printf("closest jets %d %d dR = %f",ich,ifu,dR);
446
447	chJet->SetClosestJet(fullJet,dR);
448	fullJet->SetClosestJet(chJet,dR);
449
450	}
451	}
452	}
453	}
454
455	fMatchingDone = kTRUE;
456
457	}
458
459	//________________________________________________________________________
460	void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndex() {
461
462	// take each full jet and set the index of the charged jet with the largest shared charged fraction
463
464	const Int_t nJetsCh = GetNJets(fContainerCharged);
465	const Int_t nJetsFull = GetNJets(fContainerFull);
466	faFullFracIndex.Set(nJetsFull+1);
467	faFullFracIndex.Reset(-1);
468
469	AliJetContainer *cont = GetJetContainer(fContainerFull);
470	Float_t radius = cont->GetJetRadius();
471
472	for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
473	Double_t frac = 0.;
474	Double_t dist = 10.;
475	AliEmcalJet jetFull = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ifu, fContainerFull));
476	if(!jetFull) {
477	faFullFracIndex.SetAt(-1,ifu);
478	continue;
479	}
480
481	for(Int_t ich = 0; ich<nJetsCh; ich++) {
482	AliEmcalJet jetCh = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ich, fContainerCharged));
483	if(!jetCh)
484	continue;
485	Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
486	dist = GetDeltaR(jetFull,jetCh);
487	if(tmpFrac>frac && dist<radius) {
488	frac = tmpFrac;
489	faFullFracIndex.SetAt(ich,ifu);
490	}
491	}
492
493	}
494
495	}
496
497	//________________________________________________________________________
498	void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndexMC() {
499
500	// take each full jet and set the index of the charged jet with the largest shared charged fraction
501
502	const Int_t nJetsCh = GetNJets(fContainerChargedMC);
503	const Int_t nJetsFull = GetNJets(fContainerFullMC);
504	faFullFracIndexMC.Set(nJetsFull);
505	faFullFracIndexMC.Reset(-1);
506
507	AliJetContainer *cont = GetJetContainer(fContainerFullMC);
508	Float_t radius = cont->GetJetRadius();
509
510	for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
511	Double_t frac = 0.;
512	Double_t dist = 10.;
513	AliEmcalJet jetFull = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ifu, fContainerFullMC));
514	if(!jetFull) {
515	faFullFracIndexMC.SetAt(-1,ifu);
516	continue;
517	}
518
519	for(Int_t ich = 0; ich<nJetsCh; ich++) {
520	AliEmcalJet jetCh = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ich, fContainerChargedMC));
521	if(!jetCh)
522	continue;
523	Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
524	dist = GetDeltaR(jetFull,jetCh);
525	if(tmpFrac>frac && dist<radius) {
526	frac = tmpFrac;
527	faFullFracIndexMC.SetAt(ich,ifu);
528	}
529	}
530
531	}
532
533	}
534
6ab30d5f	535	//_______________________________________________________________________
	536	AliEmcalJet* AliAnalysisTaskEmcalDiJetBase::GetLeadingJetOppositeHemisphere(const Int_t type, const Int_t typea, const AliEmcalJet *jetTrig) {
	537
	538	// Get leading jet in opposite hemisphere from trigger jet
	539	// type = correlation type
	540	// typea = container of associated jets
	541
	542	Int_t nJetsAssoc = GetNJets(typea);
	543	Double_t ptLead = -999;
	544	Int_t iJetLead = -1;
	545	for(Int_t ija=0; ija<nJetsAssoc; ija++) {
	546
	547	AliEmcalJet *jetAssoc = NULL;
	548	if(type==0) {
	549	jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
	550	if(TMath::Abs(jetAssoc->Eta())>0.5)
	551	jetAssoc = NULL;
	552	}
	553	else
	554	jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
	555
	556	if(!jetAssoc)
	557	continue;
	558
	559	Double_t dPhi = GetDeltaPhi(jetTrig,jetAssoc);
	560	Double_t phiMin = 0.5*TMath::Pi();
	561	Double_t phiMax = 1.5*TMath::Pi();
3c709670	562	if(dPhi<phiMin \|\| dPhi>phiMax)
6ab30d5f	563	continue;
	564
	565	Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
	566
	567	if(jetAssocPt>ptLead) {
	568	ptLead = jetAssocPt;
	569	iJetLead = ija;
	570	}
	571
	572	}
	573
8e49a788	574	AliEmcalJet *jetAssocLead = NULL;
	575	if(iJetLead>-1)
	576	jetAssocLead = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead, typea));
6ab30d5f	577
	578	return jetAssocLead;
	579
	580	}
	581
8e49a788	582	//_______________________________________________________________________
8e49a788	583	AliEmcalJet* AliAnalysisTaskEmcalDiJetBase::GetSecondLeadingJetOppositeHemisphere(const Int_t type, const Int_t typea, const AliEmcalJet *jetTrig) {
6e8b6371	584
8e49a788	585	// Get leading jet in opposite hemisphere from trigger jet
	586	// type = correlation type
	587	// typea = container of associated jets
6e8b6371	588
8e49a788	589	Int_t nJetsAssoc = GetNJets(typea);
	590	Double_t ptLead = -999;
	591	Int_t iJetLead = -1;
	592	Int_t iJetLead2 = -1;
	593	for(Int_t ija=0; ija<nJetsAssoc; ija++) {
	594
	595	AliEmcalJet *jetAssoc = NULL;
	596	if(type==0) {
	597	jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
	598	if(TMath::Abs(jetAssoc->Eta())>0.5)
	599	jetAssoc = NULL;
6e8b6371	600	}
8e49a788	601	else
8e49a788	602	jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
6e8b6371	603
8e49a788	604	if(!jetAssoc)
8e49a788	605	continue;
6e8b6371	606
8e49a788	607	Double_t dPhi = GetDeltaPhi(jetTrig,jetAssoc);
	608	Double_t phiMin = 0.5*TMath::Pi();
	609	Double_t phiMax = 1.5*TMath::Pi();
	610	if(dPhi<phiMin \|\| dPhi>phiMax)
	611	continue;
6e8b6371	612
8e49a788	613	Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
6e8b6371	614
8e49a788	615	if(jetAssocPt>ptLead) {
	616	iJetLead2 = iJetLead;
	617	ptLead = jetAssocPt;
	618	iJetLead = ija;
	619	}
6e8b6371	620
6e8b6371	621	}
6e8b6371	622
8e49a788	623	AliEmcalJet *jetAssocLead2 = NULL;
	624	if(iJetLead2>-1)
	625	jetAssocLead2 = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead2, typea));
6e8b6371	626
8e49a788	627	return jetAssocLead2;
6e8b6371	628
6e8b6371	629	}
	630
	631	//________________________________________________________________________
	632	Bool_t AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects() {
	633	//
8e49a788	634	// Retrieve objects from event.
6e8b6371	635	//
6e8b6371	636
9239b066	637	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
6e8b6371	638	return kFALSE;
	639
	640	if(fRhoType==0) {
	641	fRhoFullVal = 0.;
	642	fRhoChVal = 0.;
	643	}
	644	if(fRhoType==1) {
	645	fRhoFullVal = GetRhoVal(fContainerFull);
	646	fRhoChVal = GetRhoVal(fContainerCharged);
	647	}
	648
6e8b6371	649	return kTRUE;
6e8b6371	650	}
	651
	652	//_______________________________________________________________________
	653	void AliAnalysisTaskEmcalDiJetBase::Terminate(Option_t *)
	654	{
	655	// Called once at the end of the analysis.
	656	}