ana update salvatore

[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliAnalysisTaskSAJF.cxx

// $Id$
//
// Jet analysis task (S.Aiola).
//
//

#include <TObject.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TRandom3.h>
#include <TParameter.h>

#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliVEventHandler.h"
#include "AliLog.h"

#include "AliAnalysisTaskSAJF.h"

ClassImp(AliAnalysisTaskSAJF)

//________________________________________________________________________
AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() : 
  AliAnalysisTaskSE("AliAnalysisTaskSAJF"),
  fAnaType(kTPC),
  fMinEta(-0.9),
  fMaxEta(0.9),
  fMinPhi(-10),
  fMaxPhi(10),
  fJetRadius(0.4),
  fMinRC2LJ(1.0),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fEmbJetsName("EmbJets"),
  fRandTracksName("TracksRandomized"),
  fRandCaloName("CaloClustersRandomized"),
  fRhoName("Rho"),
  fNbins(500),
  fMinPt(0),
  fMaxPt(250),
  fPtCut(0.15),
  fPtBiasJetTrack(10),
  fPtBiasJetClus(10),
  fTracks(0),
  fCaloClusters(0),
  fJets(0),
  fEmbJets(0),
  fRandTracks(0),
  fRandCaloClusters(0),
  fCent(0),
  fCentBin(-1),
  fRho(0),
  fOutput(0),
  fHistCentrality(0),
  fHistJetPhiEta(0),
  fHistRhoVSleadJetPt(0),
  fHistRCPhiEta(0),
  fHistRCPtExLJVSDPhiLJ(0),
  fHistRhoVSRCPt(0),
  fHistEmbJetPhiEta(0),
  fHistEmbPartPhiEta(0),
  fHistRhoVSEmbBkg(0)

{
  // Default constructor.

  for (Int_t i = 0; i < 4; i++) {
    fHistJetsPt[i] = 0;
    fHistJetsPtTrack[i] = 0;
    fHistJetsPtClus[i] = 0;
    fHistJetsPtNonBias[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHist2LeadingJetPt[i] = 0;
    fHistJetsNEFvsPt[i] = 0;
    fHistJetsZvsPt[i] = 0;
    fHistTracksPtLJ[i] = 0;
    fHistClusEtLJ[i] = 0;
    fHistTracksPtBkg[i] = 0;
    fHistClusEtBkg[i] = 0;
    fHistRho[i] = 0;
    fHistCorrJetsPt[i] = 0;
    fHistCorrJetsPtClus[i] = 0;
    fHistCorrJetsPtTrack[i] = 0;
    fHistCorrJetsPtNonBias[i] = 0;
    fHistCorrLeadingJetPt[i] = 0;
    fHistRCPt[i] = 0;
    fHistRCPtExLJ[i] = 0;
    fHistRCPtRand[i] = 0;
    fHistDeltaPtRC[i] = 0;
    fHistDeltaPtRCExLJ[i] = 0;
    fHistDeltaPtRCRand[i] = 0;
    fHistEmbJets[i] = 0;
    fHistEmbPart[i] = 0;
    fHistDeltaPtEmb[i] = 0;
  }

  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class()); 
}

//________________________________________________________________________
AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) : 
  AliAnalysisTaskSE(name),
  fAnaType(kTPC),
  fMinEta(-0.9),
  fMaxEta(0.9),
  fMinPhi(-10),
  fMaxPhi(10),
  fJetRadius(0.4),
  fMinRC2LJ(1.0),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fJetsName("Jets"),
  fEmbJetsName("EmbJets"),
  fRandTracksName("TracksRandomized"),
  fRandCaloName("CaloClustersRandomized"),
  fRhoName("Rho"),
  fNbins(500),
  fMinPt(0),
  fMaxPt(250),
  fPtCut(0.15),
  fPtBiasJetTrack(10),
  fPtBiasJetClus(10),
  fTracks(0),
  fCaloClusters(0),
  fJets(0),
  fEmbJets(0),
  fRandTracks(0),
  fRandCaloClusters(0),
  fCent(0),
  fCentBin(-1),
  fRho(0),
  fOutput(0),
  fHistCentrality(0),
  fHistJetPhiEta(0),
  fHistRhoVSleadJetPt(0),
  fHistRCPhiEta(0),
  fHistRCPtExLJVSDPhiLJ(0),
  fHistRhoVSRCPt(0),
  fHistEmbJetPhiEta(0),
  fHistEmbPartPhiEta(0),
  fHistRhoVSEmbBkg(0)
{
  // Standard constructor.

  for (Int_t i = 0; i < 4; i++) {
    fHistJetsPt[i] = 0;
    fHistJetsPtTrack[i] = 0;
    fHistJetsPtClus[i] = 0;
    fHistJetsPtNonBias[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHist2LeadingJetPt[i] = 0;
    fHistJetsNEFvsPt[i] = 0;
    fHistJetsZvsPt[i] = 0;
    fHistTracksPtLJ[i] = 0;
    fHistClusEtLJ[i] = 0;
    fHistTracksPtBkg[i] = 0;
    fHistClusEtBkg[i] = 0;
    fHistRho[i] = 0;
    fHistCorrJetsPt[i] = 0;
    fHistCorrJetsPtClus[i] = 0;
    fHistCorrJetsPtTrack[i] = 0;
    fHistCorrJetsPtNonBias[i] = 0;
    fHistCorrLeadingJetPt[i] = 0;
    fHistRCPt[i] = 0;
    fHistRCPtExLJ[i] = 0;
    fHistRCPtRand[i] = 0;
    fHistDeltaPtRC[i] = 0;
    fHistDeltaPtRCExLJ[i] = 0;
    fHistDeltaPtRCRand[i] = 0;
    fHistEmbJets[i] = 0;
    fHistEmbPart[i] = 0;
    fHistDeltaPtEmb[i] = 0;
  }

  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class()); 
}

//________________________________________________________________________
AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
{
  // Destructor
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::UserCreateOutputObjects()
{
  // Create histograms

  const Float_t binWidth = (fMaxPt - fMinPt) / fNbins;
  
  OpenFile(1);
  fOutput = new TList();
  fOutput->SetOwner(); 
  
  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", fNbins, 0, 100);
  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
  fHistCentrality->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistCentrality);

  fHistJetPhiEta = new TH2F("fHistJetPhiEta","Phi-Eta distribution of jets", 20, -2, 2, 32, 0, 6.4);
  fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
  fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistJetPhiEta);

  fHistRhoVSleadJetPt = new TH2F("fHistRhoVSleadJetPt","fHistRhoVSleadJetPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
  fHistRhoVSleadJetPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
  fHistRhoVSleadJetPt->GetYaxis()->SetTitle("Leading jet p_{T} [GeV/c]");
  fOutput->Add(fHistRhoVSleadJetPt);

  fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 20, -2, 2, 32, 0, 6.4);
  fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
  fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistRCPhiEta);

  fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinPt, fMaxPt, 128, -1.6, 4.8);
  fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
  fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
  fOutput->Add(fHistRCPtExLJVSDPhiLJ);

  fHistRhoVSRCPt = new TH2F("fHistRhoVSRCPt","fHistRhoVSRCPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
  fHistRhoVSRCPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
  fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
  fOutput->Add(fHistRhoVSRCPt);

  fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 20, -2, 2, 32, 0, 6.4);
  fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
  fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistEmbJetPhiEta);

  fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 20, -2, 2, 32, 0, 6.4);
  fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
  fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistEmbPartPhiEta);

  fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
  fHistRhoVSEmbBkg->GetXaxis()->SetTitle("rho * area [GeV/c]");
  fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
  fOutput->Add(fHistRhoVSEmbBkg);

  TString histname;

  for (Int_t i = 0; i < 4; i++) {
    histname = "fHistJetsPt_";
    histname += i;
    fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
    fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistJetsPt[i]);

    if (fAnaType == kEMCAL) {
      histname = "fHistJetsPtClus_";
      histname += i;
      fHistJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
      fHistJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
      fHistJetsPtClus[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistJetsPtClus[i]);
    }

    histname = "fHistJetsPtTrack_";
    histname += i;
    fHistJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
    fHistJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistJetsPtTrack[i]);

    histname = "fHistJetsPtNonBias_";
    histname += i;
    fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
    fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistJetsPtNonBias[i]);

    histname = "fHistLeadingJetPt_";
    histname += i;
    fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
    fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
    fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistLeadingJetPt[i]);

    histname = "fHist2LeadingJetPt_";
    histname += i;
    fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
    fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
    fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHist2LeadingJetPt[i]);

    histname = "fHistJetsZvsPt_";
    histname += i;
    fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
    fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
    fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
    fOutput->Add(fHistJetsZvsPt[i]);

    if (fAnaType == kEMCAL) {
      histname = "fHistJetsNEFvsPt_";
      histname += i;
      fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
      fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
      fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
      fOutput->Add(fHistJetsNEFvsPt[i]);

      histname = "fHistClusEtLJ_";
      histname += i;
      fHistClusEtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
      fHistClusEtLJ[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
      fHistClusEtLJ[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistClusEtLJ[i]);
      
      histname = "fHistClusEtBkg_";
      histname += i;
      fHistClusEtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
      fHistClusEtBkg[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
      fHistClusEtBkg[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistClusEtBkg[i]);
    }

    histname = "fHistTracksPtLJ_";
    histname += i;
    fHistTracksPtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
    fHistTracksPtLJ[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistTracksPtLJ[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistTracksPtLJ[i]);

    histname = "fHistTracksPtBkg_";
    histname += i;
    fHistTracksPtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
    fHistTracksPtBkg[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistTracksPtBkg[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistTracksPtBkg[i]);

    histname = "fHistRho_";
    histname += i;
    fHistRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt * 2);
    fHistRho[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistRho[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistRho[i]);

    histname = "fHistCorrJetsPt_";
    histname += i;
    fHistCorrJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
    fHistCorrJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistCorrJetsPt[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistCorrJetsPt[i]);

    if (fAnaType == kEMCAL) {
      histname = "fHistCorrJetsPtClus_";
      histname += i;
      fHistCorrJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
      fHistCorrJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
      fHistCorrJetsPtClus[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistCorrJetsPtClus[i]);
    }

    histname = "fHistCorrJetsPtTrack_";
    histname += i;
    fHistCorrJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
    fHistCorrJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistCorrJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistCorrJetsPtTrack[i]);

    histname = "fHistCorrJetsPtNonBias_";
    histname += i;
    fHistCorrJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
    fHistCorrJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistCorrJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistCorrJetsPtNonBias[i]);

    histname = "fHistCorrLeadingJetPt_";
    histname += i;
    fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
    fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistCorrLeadingJetPt[i]);
    
    histname = "fHistRCPt_";
    histname += i;
    fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt * 2);
    fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
    fHistRCPt[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistRCPt[i]);

    histname = "fHistRCPtExLJ_";
    histname += i;
    fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt * 2);
    fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
    fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistRCPtExLJ[i]);

    histname = "fHistRCPtRand_";
    histname += i;
    fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt * 2);
    fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
    fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistRCPtRand[i]);

    histname = "fHistDeltaPtRC_";
    histname += i;
    fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
    fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
    fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistDeltaPtRC[i]);

    histname = "fHistDeltaPtRCExLJ_";
    histname += i;
    fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
    fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
    fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistDeltaPtRCExLJ[i]);

    histname = "fHistDeltaPtRCRand_";
    histname += i;
    fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
    fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
    fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistDeltaPtRCRand[i]);

    histname = "fHistEmbJets_";
    histname += i;
    fHistEmbJets[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
    fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet p_{T} [GeV/c]");
    fHistEmbJets[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistEmbJets[i]);

    histname = "fHistEmbPart_";
    histname += i;
    fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
    fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle p_{T} [GeV/c]");
    fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistEmbPart[i]);

    histname = "fHistDeltaPtEmb_";
    histname += i;
    fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
    fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
    fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistDeltaPtEmb[i]);
  }

  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::RetrieveEventObjects()
{
  if (strcmp(fCaloName,"") && fAnaType == kEMCAL) {
    fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
    if (!fCaloClusters) {
      AliWarning(Form("Could not retrieve clusters %s!", fCaloName.Data())); 
    }
  }

  if (strcmp(fTracksName,"")) {
    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
    if (!fTracks) {
      AliWarning(Form("Could not retrieve tracks %s!", fTracksName.Data())); 
    }
  }

  if (strcmp(fJetsName,"")) {
    fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
    if (!fJets) {
      AliWarning(Form("Could not retrieve jets %s!", fJetsName.Data())); 
    }
  }

  if (strcmp(fEmbJetsName,"")) {
    fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
    if (!fEmbJets) {
      AliWarning(Form("Could not retrieve emb jets %s!", fEmbJetsName.Data())); 
    }
  }

  if (strcmp(fRandCaloName,"") && fAnaType == kEMCAL) {
    fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
    if (!fRandCaloClusters) {
      AliWarning(Form("Could not retrieve randomized clusters %s!", fRandCaloName.Data())); 
    }
  }

  if (strcmp(fRandTracksName,"")) {
    fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
    if (!fRandTracks) {
      AliWarning(Form("Could not retrieve randomized tracks %s!", fRandTracksName.Data())); 
    }
  }

  AliCentrality *aliCent = InputEvent()->GetCentrality();
  if (aliCent) {
    fCent = aliCent->GetCentralityPercentile("V0M");
    if      (fCent >=  0 && fCent <   10) fCentBin = 0;
    else if (fCent >= 10 && fCent <   30) fCentBin = 1;
    else if (fCent >= 30 && fCent <   50) fCentBin = 2;
    else if (fCent >= 50 && fCent <= 100) fCentBin = 3; 
    else {
      AliWarning(Form("Negative centrality: %f. Assuming 99", fCent));
      fCentBin = 3;
    }
  }
  else {
    AliWarning(Form("Could not retrieve centrality information! Assuming 99"));
    fCentBin = 3;
  }

  fRho = -1;
  
  if (strcmp(fRhoName,"")) {
    TParameter<Double_t> *rho = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhoName));
  
    if (rho) {
      fRho = rho->GetVal();
    }
    else {
      AliWarning(Form("Could not retrieve rho %s!", fRhoName.Data()));
    }
  }

  fVertex[0] = 0;
  fVertex[1] = 0;
  fVertex[2] = 0;
  InputEvent()->GetPrimaryVertex()->GetXYZ(fVertex);
}

//________________________________________________________________________
AliVTrack* AliAnalysisTaskSAJF::GetTrack(const Int_t i) const
{
  if (fTracks)
    return dynamic_cast<AliVTrack*>(fTracks->At(i));
  else
    return 0;
}

//________________________________________________________________________
Int_t AliAnalysisTaskSAJF::GetNumberOfTracks() const
{
  if (fTracks)
    return fTracks->GetEntriesFast();
  else
    return 0;
}

//________________________________________________________________________
AliVCluster* AliAnalysisTaskSAJF::GetCaloCluster(const Int_t i) const
{ 
  if (fCaloClusters)
    return dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
  else
    return 0;
}

//________________________________________________________________________
Int_t AliAnalysisTaskSAJF::GetNumberOfCaloClusters() const
{ 
  if (fCaloClusters)
    return fCaloClusters->GetEntriesFast();
  else
    return 0;
}

//________________________________________________________________________
AliEmcalJet* AliAnalysisTaskSAJF::GetJet(const Int_t i) const
{
  if (fJets)
    return dynamic_cast<AliEmcalJet*>(fJets->At(i));
  else
    return 0;
}

//________________________________________________________________________
Int_t AliAnalysisTaskSAJF::GetNumberOfJets() const
{
  if (fJets)
    return fJets->GetEntriesFast();
  else
    return 0;
}

//________________________________________________________________________
AliEmcalJet* AliAnalysisTaskSAJF::GetEmbJet(const Int_t i) const
{
  if (fEmbJets)
    return dynamic_cast<AliEmcalJet*>(fEmbJets->At(i));
  else
    return 0;
}

//________________________________________________________________________
Int_t AliAnalysisTaskSAJF::GetNumberOfEmbJets() const
{
  if (fEmbJets)
    return fEmbJets->GetEntriesFast();
  else
    return 0;
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::FillHistograms()
{
  const Float_t A = fJetRadius * fJetRadius * TMath::Pi();

  Int_t maxJetIndex  = -1;
  Int_t max2JetIndex = -1;

  // ************
  // General histograms
  // _________________________________

  DoJetLoop(maxJetIndex, max2JetIndex);
  
  if (maxJetIndex < 0) 
    return;

  AliEmcalJet* jet = GetJet(maxJetIndex);
  if (!jet) 
    return;

  fHistCentrality->Fill(fCent);

  fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
  fHistRhoVSleadJetPt->Fill(fRho * jet->Area(), jet->Pt());

  jet->SortConstituents();
  
  AliEmcalJet* jet2 = 0;
  if (max2JetIndex >= 0)
    jet2 = GetJet(max2JetIndex);

  if (jet2) {
    fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
    jet2->SortConstituents();
  }

  fHistRho[fCentBin]->Fill(fRho);

  Float_t maxJetCorrPt = jet->Pt() - fRho * jet->Area();
  if (maxJetCorrPt > 0)
    fHistCorrLeadingJetPt[fCentBin]->Fill(maxJetCorrPt);
  
  DoTrackLoop(maxJetIndex);

  if (fAnaType == kEMCAL)
    DoClusterLoop(maxJetIndex);

  // ************
  // Random cones
  // _________________________________
  
  // Simple random cones
  Float_t RCpt = 0;
  Float_t RCeta = 0;
  Float_t RCphi = 0;
  GetRigidCone(RCpt, RCeta, RCphi, kFALSE, 0);
  if (RCpt > 0) {
    fHistRCPt[fCentBin]->Fill(RCpt / A);
    fHistDeltaPtRC[fCentBin]->Fill(RCpt - A * fRho);
  }
  
  // Random cones far from leading jet
  Float_t RCptExLJ = 0;
  Float_t RCetaExLJ = 0;
  Float_t RCphiExLJ = 0;
  GetRigidCone(RCptExLJ, RCetaExLJ, RCphiExLJ, kFALSE, jet);
  if (RCptExLJ > 0) {
    fHistRCPhiEta->Fill(RCetaExLJ, RCphiExLJ);
    fHistRhoVSRCPt->Fill(fRho, RCptExLJ / A);

    Float_t dphi = RCphiExLJ - jet->Phi();
    if (dphi > 4.8) dphi -= TMath::Pi() * 2;
    if (dphi < -1.6) dphi += TMath::Pi() * 2; 
    fHistRCPtExLJVSDPhiLJ->Fill(RCptExLJ, dphi);
    
    fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / A);
    fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - A * fRho);
  }

  // Random cones with randomized particles
  Float_t RCptRand = 0;
  Float_t RCetaRand = 0;
  Float_t RCphiRand = 0;
  GetRigidCone(RCptRand, RCetaRand, RCphiRand, kTRUE, 0, fRandTracks, fRandCaloClusters);
  if (RCptRand > 0) {
    fHistRCPtRand[fCentBin]->Fill(RCptRand / A);
    fHistDeltaPtRCRand[fCentBin]->Fill(RCptRand - A * fRho);
  }

  // ************
  // Embedding
  // _________________________________

  AliEmcalJet *embJet  = 0;
  TObject     *maxPart = 0;

  DoEmbJetLoop(embJet, maxPart);

  if (embJet && maxPart) {
    Float_t maxEmbPartPt = 0;
    Float_t maxEmbPartEta = 0;
    Float_t maxEmbPartPhi = 0;

    AliVCluster *cluster = dynamic_cast<AliVCluster*>(maxPart);
    if (cluster) {
      TLorentzVector nPart;
      cluster->GetMomentum(nPart, fVertex);
      Float_t pos[3];
      cluster->GetPosition(pos);
      TVector3 clusVec(pos);
      
      maxEmbPartPt = nPart.Et();
      maxEmbPartEta = clusVec.Eta();
      maxEmbPartPhi = clusVec.Phi();
    }
    else {
      AliVTrack *track = dynamic_cast<AliVTrack*>(maxPart);
      if (track) {
	maxEmbPartPt = track->Pt();
	maxEmbPartEta = track->Eta();
	maxEmbPartPhi = track->Phi();
      }
      else {
	AliWarning(Form("%s - Embedded particle type not recognized (neither AliVCluster nor AliVTrack)!", GetName()));
	return;
      }
    }
    fHistEmbJets[fCentBin]->Fill(embJet->Pt());
    fHistEmbPart[fCentBin]->Fill(maxEmbPartPt);
    fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
    fHistEmbPartPhiEta->Fill(maxEmbPartEta, maxEmbPartPhi);

    fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRho - maxEmbPartPt);
    fHistRhoVSEmbBkg->Fill(embJet->Area() * fRho, embJet->Pt() - maxEmbPartPt);
  }
  else {
    AliWarning(Form("%s - Embedded particle not found in the event!", GetName()));
  }
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
{
  Int_t njets = GetNumberOfJets();

  Float_t maxJetPt = 0;
  Float_t max2JetPt = 0;
  for (Int_t ij = 0; ij < njets; ij++) {

    AliEmcalJet* jet = GetJet(ij);

    if (!jet) {
      AliError(Form("Could not receive jet %d", ij));
      continue;
    }  
    
    if (jet->Pt() <= 0)
      continue;

    if (!AcceptJet(jet))
      continue;

    Float_t corrPt = jet->Pt() - fRho * jet->Area();

    fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
    fHistCorrJetsPtNonBias[fCentBin]->Fill(corrPt);

    if (jet->MaxTrackPt() > fPtBiasJetTrack) {
      fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
      fHistCorrJetsPtTrack[fCentBin]->Fill(corrPt);
    }
    
    if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus) {
      fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
      fHistCorrJetsPtClus[fCentBin]->Fill(corrPt);
    }
    
    if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet->MaxClusterPt() < fPtBiasJetClus))
	continue;

    fHistJetsPt[fCentBin]->Fill(jet->Pt());
    fHistCorrJetsPt[fCentBin]->Fill(corrPt);

    fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());

    if (fAnaType == kEMCAL)
      fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt());

    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
      AliVTrack *track = jet->TrackAt(it, fTracks);
      if (track)
	fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt());
    }

    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
      AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);

      if (cluster) {
	TLorentzVector nPart;
	cluster->GetMomentum(nPart, fVertex);
	fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), corrPt);
      }
    }

    if (maxJetPt < corrPt) {
      max2JetPt = maxJetPt;
      max2JetIndex = maxJetIndex;
      maxJetPt = corrPt;
      maxJetIndex = ij;
    }
    else if (max2JetPt < corrPt) {
      max2JetPt = corrPt;
      max2JetIndex = ij;
    }
  } //jet loop 
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &maxPart)
{
  TLorentzVector *maxClusVect = new TLorentzVector();

  Int_t nembjets = GetNumberOfEmbJets();

  for (Int_t ij = 0; ij < nembjets; ij++) {
      
    AliEmcalJet* jet = GetEmbJet(ij);
      
    if (!jet) {
      AliError(Form("Could not receive jet %d", ij));
      continue;
    } 
      
    if (jet->Pt() <= 0)
	continue;
 
    //if (!AcceptJet(jet))
    //continue;

    AliVTrack *maxTrack = 0;

    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
      AliVTrack *track = jet->TrackAt(it, fTracks);
      
      if (!track) continue;
      
      if (!maxTrack || track->Pt() > maxTrack->Pt())
	maxTrack = track;
    }
    
    AliVCluster *maxClus = 0;

    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
      AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
      
      if (!cluster) continue;
      
      TLorentzVector nPart;
      cluster->GetMomentum(nPart, fVertex);
      
      if (!maxClus || nPart.Et() > maxClusVect->Et()) {
	new (maxClusVect) TLorentzVector(nPart);
	maxClus = cluster;
      } 
    }

    if ((maxClus && maxTrack && maxClusVect->Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
      if (maxClus->Chi2() == 100) {
	maxPart = maxClus;
	embJet = jet;
	delete maxClusVect;
      	return;
      }
    }
    else if (maxTrack) {
      if (maxTrack->GetLabel() == 100) {
	maxPart = maxTrack;
	embJet = jet;
	delete maxClusVect;
	return;
      }
    }
  }

  delete maxClusVect;
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::DoTrackLoop(Int_t maxJetIndex)
{ 
  AliEmcalJet* jet = 0;
  if (maxJetIndex >= 0)
    jet = GetJet(maxJetIndex);

  Int_t ntracks = GetNumberOfTracks();

  for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
    AliVTrack* track = GetTrack(iTracks);         
    if(!track) {
      AliError(Form("Could not retrieve track %d",iTracks)); 
      continue; 
    }
    
    if (!AcceptTrack(track)) continue;

    if (jet && IsJetTrack(jet, iTracks)) {
      fHistTracksPtLJ[fCentBin]->Fill(track->Pt());
    }
    else {
      fHistTracksPtBkg[fCentBin]->Fill(track->Pt());
    }
  }
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::DoClusterLoop(Int_t maxJetIndex)
{
  AliEmcalJet* jet = 0;
  if (maxJetIndex >= 0)
    jet = GetJet(maxJetIndex);

  Int_t nclusters =  GetNumberOfCaloClusters();
  for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
    AliVCluster* cluster = GetCaloCluster(iClusters);
    if (!cluster) {
      AliError(Form("Could not receive cluster %d", iClusters));
      continue;
    }  
    
    if (!(cluster->IsEMCAL())) continue;
    
    if (!AcceptCluster(cluster)) continue;

    TLorentzVector nPart;
    cluster->GetMomentum(nPart, fVertex);

    if (jet && IsJetCluster(jet, iClusters)) {
      fHistClusEtLJ[fCentBin]->Fill(nPart.Et());
    }
    else {
      fHistClusEtBkg[fCentBin]->Fill(nPart.Et());
    }
  } //cluster loop 
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::Init()
{
  if (fAnaType == kTPC) {
    fMinEta = -0.9;
    fMaxEta = 0.9;
    fMinPhi = -10;
    fMaxPhi = 10;
  }
  else if (fAnaType == kEMCAL) {
    fMinEta = -0.7;
    fMaxEta = 0.7;
    fMinPhi = 80 * TMath::DegToRad();
    fMaxPhi = 180 * TMath::DegToRad();
  }
  else {
    AliWarning("Analysis type not recognized! Assuming kTPC...");
    fAnaType = kTPC;
    Init();
  }

  const Float_t semiDiag = TMath::Sqrt((fMaxPhi - fMinPhi) * (fMaxPhi - fMinPhi) + (fMaxEta - fMinEta) * (fMaxEta - fMinEta)) / 2;
  if (fMinRC2LJ > semiDiag * 0.5) {
    AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. Will use fMinRC2LJ = %f", fMinRC2LJ, semiDiag * 0.5));
    fMinRC2LJ = semiDiag * 0.5;
  }
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::GetRigidCone(Float_t &pt, Float_t &eta, Float_t &phi, Bool_t acceptMC,
				       AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
{
  if (!tracks)
    tracks = fTracks;

  if (!clusters)
    clusters = fCaloClusters;

  eta = 0;
  phi = 0;
  pt = 0;

  Float_t LJeta = 999;
  Float_t LJphi = 999;

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

  Float_t maxEta = fMaxEta - fJetRadius;
  Float_t minEta = fMinEta + fJetRadius;
  Float_t maxPhi = fMaxPhi - fJetRadius;
  Float_t minPhi = fMinPhi + fJetRadius;

  if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
  if (minPhi < 0) minPhi = 0;
  
  Float_t dLJ = 0;
  Int_t repeats = 0;
  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);

  if (repeats == 999)
    return;

  if (fAnaType == kEMCAL) {
    Int_t nclusters =  clusters->GetEntriesFast();
    for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
      AliVCluster* cluster = dynamic_cast<AliVCluster*>(clusters->At(iClusters));
      if (!cluster) {
	AliError(Form("Could not receive cluster %d", iClusters));
	continue;
      }  
      
      if (!(cluster->IsEMCAL())) continue;
      
      if (!AcceptCluster(cluster, acceptMC)) continue;
      
      TLorentzVector nPart;
      cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
      
      Float_t pos[3];
      cluster->GetPosition(pos);
      TVector3 clusVec(pos);
      
      Float_t d = TMath::Sqrt((clusVec.Eta() - eta) * (clusVec.Eta() - eta) + (clusVec.Phi() - phi) * (clusVec.Phi() - phi));
      
      if (d <= fJetRadius)
	pt += nPart.Pt();
    }
  }

  Int_t ntracks = tracks->GetEntriesFast();
  for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
    AliVTrack* track = dynamic_cast<AliVTrack*>(tracks->At(iTracks));         
    if(!track) {
      AliError(Form("Could not retrieve track %d",iTracks)); 
      continue; 
    }
    
    if (!AcceptTrack(track, acceptMC)) continue;

    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 <= fJetRadius)
      pt += track->Pt();
  }
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSAJF::IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted) const
{
  for (Int_t i = 0; i < jet->GetNumberOfTracks(); i++) {
    Int_t ijettrack = jet->TrackAt(i);
    if (sorted && ijettrack > itrack)
      return kFALSE;
    if (ijettrack == itrack)
      return kTRUE;
  }
  return kFALSE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSAJF::IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted) const
{
  for (Int_t i = 0; i < jet->GetNumberOfClusters(); i++) {
    Int_t ijetclus = jet->ClusterAt(i);
    if (sorted && ijetclus > iclus)
      return kFALSE;
    if (ijetclus == iclus)
      return kTRUE;
  }
  return kFALSE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSAJF::AcceptJet(Float_t eta, Float_t phi) const
{
  return (Bool_t)(eta > fMinEta + fJetRadius && eta < fMaxEta - fJetRadius && phi > fMinPhi + fJetRadius && phi < fMaxPhi - fJetRadius);
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSAJF::AcceptJet(AliEmcalJet *jet) const
{
  return AcceptJet(jet->Eta(), jet->Phi());
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSAJF::AcceptCluster(AliVCluster* clus, Bool_t acceptMC) const
{
  if (!acceptMC && clus->Chi2() == 100)
    return kFALSE;

  TLorentzVector nPart;
  clus->GetMomentum(nPart, const_cast<Double_t*>(fVertex));

  if (nPart.Et() < fPtCut)
    return kFALSE;

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSAJF::AcceptTrack(AliVTrack* track, Bool_t acceptMC) const
{
  if (!acceptMC && track->GetLabel() == 100)
    return kFALSE;

  if (track->Pt() < fPtCut)
    return kFALSE;
  
  return (Bool_t)(track->Eta() > fMinEta && track->Eta() < fMaxEta && track->Phi() > fMinPhi && track->Phi() < fMaxPhi);
}

//________________________________________________________________________
void AliAnalysisTaskSAJF::UserExec(Option_t *) 
{
  Init();

  RetrieveEventObjects();

  if (fRho < 0) {
    AliWarning(Form("Could not retrieve rho information! Event skipped!"));
    return;
  }

  FillHistograms();
    
  // information for this iteration of the UserExec in the container
  PostData(1, fOutput);
}

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