//
// Jet mass background 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 <TRandom3.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 "AliClusterContainer.h"
#include "AliParticleContainer.h"

#include "AliAnalysisTaskEmcalJetMassBkg.h"

ClassImp(AliAnalysisTaskEmcalJetMassBkg)

//________________________________________________________________________
AliAnalysisTaskEmcalJetMassBkg::AliAnalysisTaskEmcalJetMassBkg() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMassBkg", kTRUE),
  fContainerBase(0),
  fMinRC2LJ(-1),
  fRCperEvent(10),
  fConeRadius(0.2),
  fConeMinEta(-0.9),
  fConeMaxEta(0.9),
  fConeMinPhi(0),
  fConeMaxPhi(TMath::Pi()*2),
  fJetsCont(0),
  fTracksCont(0),
  fCaloClustersCont(0),
  fh2PtVsMassRC(0),
  fpPtVsMassRC(0),
  fh2PtVsMassRCExLJDPhi(0),
  fpPtVsMassRCExLJ(0),
  fh2PtVsMassPerpConeLJ(0),
  fpPtVsMassPerpConeLJ(0),
  fh2PtVsMassPerpConeTJ(0),
  fpPtVsMassPerpConeTJ(0),
  fh2CentVsMassRC(0),
  fh2CentVsMassRCExLJ(0),
  fh2CentVsMassPerpConeLJ(0),
  fh2CentVsMassPerpConeTJ(0),
  fh2MultVsMassRC(0),
  fh2MultVsMassRCExLJ(0),
  fh2MultVsMassPerpConeLJ(0),
  fh2MultVsMassPerpConeTJ(0)
{
  // Default constructor.

  fh2PtVsMassRC            = new TH2F*[fNcentBins];
  fpPtVsMassRC             = new TProfile*[fNcentBins];
  fh2PtVsMassRCExLJDPhi    = new TH3F*[fNcentBins];
  fpPtVsMassRCExLJ         = new TProfile*[fNcentBins];
  fh2PtVsMassPerpConeLJ    = new TH2F*[fNcentBins];
  fpPtVsMassPerpConeLJ     = new TProfile*[fNcentBins];
  fh2PtVsMassPerpConeTJ    = new TH2F*[fNcentBins];
  fpPtVsMassPerpConeTJ     = new TProfile*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh2PtVsMassRC[i]             = 0;
    fpPtVsMassRC[i]              = 0;
    fh2PtVsMassRCExLJDPhi[i]     = 0;
    fpPtVsMassRCExLJ[i]          = 0;
    fh2PtVsMassPerpConeLJ[i]     = 0;
    fpPtVsMassPerpConeLJ[i]      = 0;
    fh2PtVsMassPerpConeTJ[i]     = 0;
    fpPtVsMassPerpConeTJ[i]      = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
  
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetMassBkg::AliAnalysisTaskEmcalJetMassBkg(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),  
  fContainerBase(0),
  fMinRC2LJ(-1),
  fRCperEvent(10),
  fConeRadius(0.2),
  fConeMinEta(-0.9),
  fConeMaxEta(0.9),
  fConeMinPhi(0),
  fConeMaxPhi(TMath::Pi()*2),
  fJetsCont(0),
  fTracksCont(0),
  fCaloClustersCont(0),
  fh2PtVsMassRC(0),
  fpPtVsMassRC(0),
  fh2PtVsMassRCExLJDPhi(0),
  fpPtVsMassRCExLJ(0),
  fh2PtVsMassPerpConeLJ(0),
  fpPtVsMassPerpConeLJ(0),
  fh2PtVsMassPerpConeTJ(0),
  fpPtVsMassPerpConeTJ(0),
  fh2CentVsMassRC(0),
  fh2CentVsMassRCExLJ(0),
  fh2CentVsMassPerpConeLJ(0),
  fh2CentVsMassPerpConeTJ(0),
  fh2MultVsMassRC(0),
  fh2MultVsMassRCExLJ(0),
  fh2MultVsMassPerpConeLJ(0),
  fh2MultVsMassPerpConeTJ(0)
{
  // Standard constructor.

  fh2PtVsMassRC            = new TH2F*[fNcentBins];
  fpPtVsMassRC             = new TProfile*[fNcentBins];
  fh2PtVsMassRCExLJDPhi    = new TH3F*[fNcentBins];
  fpPtVsMassRCExLJ         = new TProfile*[fNcentBins];
  fh2PtVsMassPerpConeLJ    = new TH2F*[fNcentBins];
  fpPtVsMassPerpConeLJ     = new TProfile*[fNcentBins];
  fh2PtVsMassPerpConeTJ    = new TH2F*[fNcentBins];
  fpPtVsMassPerpConeTJ     = new TProfile*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh2PtVsMassRC[i]              = 0;
    fpPtVsMassRC[i]               = 0;
    fh2PtVsMassRCExLJDPhi[i]      = 0;
    fpPtVsMassRCExLJ[i]           = 0;
    fh2PtVsMassPerpConeLJ[i]       = 0;
    fpPtVsMassPerpConeLJ[i]        = 0;
    fh2PtVsMassPerpConeTJ[i]       = 0;
    fpPtVsMassPerpConeTJ[i]        = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

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

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

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  fJetsCont         = GetJetContainer(fContainerBase);
  fTracksCont       = fJetsCont->GetParticleContainer();
  fCaloClustersCont = fJetsCont->GetClusterContainer();

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

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

  const Int_t nBinsCent  = 100;
  const Double_t minCent = 0.;
  const Double_t maxCent = 100.;

  const Int_t nBinsMult  = 400;
  const Double_t minMult = 0.;
  const Double_t maxMult = 4000.;

  fh2CentVsMassRC = new TH2F("fh2CentVsMassRC","fh2CentVsMassRC;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2CentVsMassRC);  

  fh2CentVsMassRCExLJ = new TH2F("fh2CentVsMassRCExLJ","fh2CentVsMassRCExLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2CentVsMassRCExLJ);  

  fh2CentVsMassPerpConeLJ = new TH2F("fh2CentVsMassPerpConeLJ","fh2CentVsMassPerpConeLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2CentVsMassPerpConeLJ);  

  fh2CentVsMassPerpConeTJ = new TH2F("fh2CentVsMassPerpConeTJ","fh2CentVsMassPerpConeTJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2CentVsMassPerpConeTJ); 

  fh2MultVsMassRC = new TH2F("fh2MultVsMassRC","fh2MultVsMassRC;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2MultVsMassRC);  

  fh2MultVsMassRCExLJ = new TH2F("fh2MultVsMassRCExLJ","fh2MultVsMassRCExLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2MultVsMassRCExLJ);  

  fh2MultVsMassPerpConeLJ = new TH2F("fh2MultVsMassPerpConeLJ","fh2MultVsMassPerpConeLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2MultVsMassPerpConeLJ);  

  fh2MultVsMassPerpConeTJ = new TH2F("fh2MultVsMassPerpConeTJ","fh2MultVsMassPerpConeTJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2MultVsMassPerpConeTJ); 

  TString histName = "";
  TString histTitle = "";
  for (Int_t i = 0; i < fNcentBins; i++) {
    histName = TString::Format("fh2PtVsMassRC_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
    fh2PtVsMassRC[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
    fOutput->Add(fh2PtVsMassRC[i]);

    histName = TString::Format("fpPtVsMassRC_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
    fpPtVsMassRC[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
    fOutput->Add(fpPtVsMassRC[i]);

    histName = TString::Format("fh2PtVsMassRCExLJDPhi_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
    fh2PtVsMassRCExLJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,72,-0.5*TMath::Pi(),1.5*TMath::Pi());
    fOutput->Add(fh2PtVsMassRCExLJDPhi[i]);

    histName = TString::Format("fpPtVsMassRCExLJ_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
    fpPtVsMassRCExLJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
    fOutput->Add(fpPtVsMassRCExLJ[i]);

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

    histName = TString::Format("fpPtVsMassPerpConeLJ_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
    fpPtVsMassPerpConeLJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
    fOutput->Add(fpPtVsMassPerpConeLJ[i]);

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

    histName = TString::Format("fpPtVsMassPerpConeTJ_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
    fpPtVsMassPerpConeTJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
    fOutput->Add(fpPtVsMassPerpConeTJ[i]);
  }

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

  TH1::AddDirectory(oldStatus);

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

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

  return kTRUE;
}

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

  const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
  Double_t rho = GetRhoVal(fContainerBase);
  Int_t trackMult = fTracksCont->GetNAcceptedParticles();

  TLorentzVector lvRC(0.,0.,0.,0.);
  Float_t RCpt = 0;
  Float_t RCeta = 0;
  Float_t RCphi = 0;
  Float_t RCmass = 0.;  
  for (Int_t i = 0; i < fRCperEvent; i++) {
    // Simple random cones
    lvRC.SetPxPyPzE(0.,0.,0.,0.);
    RCpt = 0;
    RCeta = 0;
    RCphi = 0;
    GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
    RCmass = lvRC.M();
    if (RCpt > 0) {
      fh2PtVsMassRC[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
      fpPtVsMassRC[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
      fh2CentVsMassRC->Fill(fCent,RCmass);
      fh2MultVsMassRC->Fill(trackMult,RCmass);
    }

    if (fJetsCont) {

      // Random cones far from leading jet(s)
      AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
      lvRC.SetPxPyPzE(0.,0.,0.,0.);
      RCpt = 0;
      RCeta = 0;
      RCphi = 0;
      GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
      RCmass = lvRC.M();
      if (RCpt > 0 && jet) {
	Float_t dphi = RCphi - jet->Phi();
	if (dphi > 1.5*TMath::Pi()) dphi -= TMath::Pi() * 2;
	if (dphi < -0.5*TMath::Pi()) dphi += TMath::Pi() * 2;
	fh2PtVsMassRCExLJDPhi[fCentBin]->Fill(RCpt - rho*rcArea,RCmass,dphi);
	fpPtVsMassRCExLJ[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
	fh2CentVsMassRCExLJ->Fill(fCent,RCmass);
	fh2MultVsMassRCExLJ->Fill(trackMult,RCmass);
      }
    }
  }//RC loop

  if(fJetsCont) {
    //cone perpendicular to leading jet
    TLorentzVector lvPC(0.,0.,0.,0.);
    Float_t PCpt = 0;
    Float_t PCeta = 0;
    Float_t PCphi = 0;
    Float_t PCmass = 0.;  
    AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
    if(jet) {
      GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
      PCmass = lvPC.M();
      if(PCpt>0.) {
	fh2PtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
	fpPtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
	fh2CentVsMassPerpConeLJ->Fill(fCent,PCmass);
	fh2MultVsMassPerpConeLJ->Fill(trackMult,PCmass);
      }
    }
    //cone perpendicular to all tagged jets
    for(int i = 0; i < fJetsCont->GetNJets();++i) {
      jet = static_cast<AliEmcalJet*>(fJetsCont->GetAcceptJet(i));
      if(!jet) continue;

      if(jet->GetTagStatus()<1)
	continue;

      lvPC.SetPxPyPzE(0.,0.,0.,0.);
      PCpt = 0;
      PCeta = 0;
      PCphi = 0;
      GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
      PCmass = lvPC.M();
      if(PCpt>0.) {
	fh2PtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
	fpPtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
	fh2CentVsMassPerpConeTJ->Fill(fCent,PCmass);
	fh2MultVsMassPerpConeTJ->Fill(trackMult,PCmass);
      }
    }//jet loop
  }

  return kTRUE;

}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetMassBkg::GetRandomCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi,
						   AliParticleContainer* tracks, AliClusterContainer* clusters,
						   AliEmcalJet *jet) const
{
  // Get rigid cone.
  lvRC.SetPxPyPzE(0.,0.,0.,0.);

  eta = -999;
  phi = -999;
  pt = 0;

  if (!tracks && !clusters)
    return;

  Float_t LJeta = 999;
  Float_t LJphi = 999;

  if (jet) {
    LJeta = jet->Eta();
    LJphi = jet->Phi();
  }

  Float_t maxEta = fConeMaxEta;
  Float_t minEta = fConeMinEta;
  Float_t maxPhi = fConeMaxPhi;
  Float_t minPhi = fConeMinPhi;

  if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
  if (minPhi < 0) minPhi = 0;
  
  Float_t dLJ = 0;
  Int_t repeats = 0;
  Bool_t reject = kTRUE;
  do {
    eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
    phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
    dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));

    repeats++;
  } while (dLJ < fMinRC2LJ && repeats < 999 && reject);

  if (repeats == 999) {
    AliWarning(Form("%s: Could not get random cone!", GetName()));
    return;
  }

  GetCone(lvRC,pt,eta,phi,tracks,clusters);


}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetMassBkg::GetCone(TLorentzVector& lvRC,Float_t &pt, Float_t eta, Float_t phi, AliParticleContainer* tracks, AliClusterContainer* clusters) const
{

  pt = 0.;
  lvRC.SetPxPyPzE(0.,0.,0.,0.);

  if (clusters) {
    AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
    while (cluster) {     
      TLorentzVector nPart;
      cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));

      Float_t cluseta = nPart.Eta();
      Float_t clusphi = nPart.Phi();
      
      if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
	clusphi += 2 * TMath::Pi();
      if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
	clusphi -= 2 * TMath::Pi();
     
      Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
      if (d <= fConeRadius) {
	pt += nPart.Pt();
	TLorentzVector lvcl(nPart.Px(),nPart.Py(),nPart.Pz(),nPart.E());
	lvRC += lvcl;
      }

      cluster = clusters->GetNextAcceptCluster();
    }
  }

  if (tracks) {
    AliVParticle* track = tracks->GetNextAcceptParticle(0); 
    while(track) { 
      Float_t tracketa = track->Eta();
      Float_t trackphi = track->Phi();
      
      if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
	trackphi += 2 * TMath::Pi();
      if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
	trackphi -= 2 * TMath::Pi();
      
      Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
      if (d <= fConeRadius) {
	pt += track->Pt();
	TLorentzVector lvtr(track->Px(),track->Py(),track->Pz(),track->E());
	lvRC += lvtr;
      }

      track = tracks->GetNextAcceptParticle(); 
    }
  }

}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetMassBkg::GetPerpCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer* tracks, AliClusterContainer* clusters, AliEmcalJet *jet) const
{
  // Get rigid cone.
  lvRC.SetPxPyPzE(0.,0.,0.,0.);

  eta = -999;
  phi = -999;
  pt = 0;

  if (!tracks && !clusters)
    return;

  if(!jet)
    return;

  Float_t LJeta = jet->Eta();
  Float_t LJphi = jet->Phi();

  eta = LJeta;
  phi = LJphi + 0.5*TMath::Pi();
  if(phi>TMath::TwoPi()) phi-=TMath::TwoPi();

  GetCone(lvRC,pt,eta,phi,tracks,clusters);
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiEMCAL()
{
  // Set default cuts for full cones

  SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
  SetConePhiLimits(1.405+fConeRadius,3.135-fConeRadius);
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiTPC()
{
  // Set default cuts for charged cones

  SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
  SetConePhiLimits(-10, 10);
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetMassBkg::ExecOnce() {

  AliAnalysisTaskEmcalJet::ExecOnce();

  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;

  if (fRCperEvent < 0) {
    Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
    Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
    fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
    if (fRCperEvent == 0)
      fRCperEvent = 1;
  }

  if (fMinRC2LJ < 0)
    fMinRC2LJ = fConeRadius * 1.5;

  const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
  if (fMinRC2LJ > maxDist) {
    AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
                    "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
    fMinRC2LJ = maxDist;
  }

}

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

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

  return kTRUE;

}

//_______________________________________________________________________
void AliAnalysisTaskEmcalJetMassBkg::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}