]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Changes from Salvatore:
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 3 Sep 2012 07:00:58 +0000 (07:00 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 3 Sep 2012 07:00:58 +0000 (07:00 +0000)
1) I splitted my jet task into a deltaPt task (AliAnalysisTaskDeltaPt) and a jet spectra task (AliAnalysisTaskSAJF)
2) I derived AliAnalysisTaskRhoBase from AliAnalysisTaskEmcalJet
3) QA rho plots are made now in AliAnalysisTaskRhoBase (so every rho implementation can benefit from them); an option can disable them as before
4) I changed my rho meth 4 (mostly for some "academic" plots for my thesis)

24 files changed:
PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJet.cxx
PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJet.h
PWGJE/EMCALJetTasks/AliAnalysisTaskRho.cxx
PWGJE/EMCALJetTasks/AliAnalysisTaskRho.h
PWGJE/EMCALJetTasks/AliAnalysisTaskRhoAverage.cxx
PWGJE/EMCALJetTasks/AliAnalysisTaskRhoAverage.h
PWGJE/EMCALJetTasks/AliAnalysisTaskRhoBase.cxx
PWGJE/EMCALJetTasks/AliAnalysisTaskRhoBase.h
PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
PWGJE/EMCALJetTasks/AliJetResponseMaker.h
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.cxx
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.h
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.h
PWGJE/EMCALJetTasks/macros/AddTaskDeltaPt.C [new file with mode: 0644]
PWGJE/EMCALJetTasks/macros/AddTaskRho.C
PWGJE/EMCALJetTasks/macros/AddTaskRhoAverage.C
PWGJE/EMCALJetTasks/macros/AddTaskRhoBase.C
PWGJE/EMCALJetTasks/macros/AddTaskSAJF.C
PWGJE/EMCALJetTasks/macros/AddTaskSAQA.C
PWGJE/PWGJEEMCALJetTasksLinkDef.h

index c464b50a3dd0634b062f9ca9b72bd7e2470c01f3..d6a59a2b02c0f446e7d74e8d29c9024c3b59e44d 100644 (file)
 
 set ( SRCS 
  EMCALJetTasks/AliAnalysisTaskEmcalJet.cxx
+ EMCALJetTasks/AliAnalysisTaskRhoBase.cxx
  EMCALJetTasks/AliAnalysisTaskRho.cxx
  EMCALJetTasks/AliAnalysisTaskRhoAverage.cxx
- EMCALJetTasks/AliAnalysisTaskRhoBase.cxx
+ EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx
  EMCALJetTasks/AliAnalysisTaskScale.cxx
  EMCALJetTasks/AliEmcalJet.cxx
  EMCALJetTasks/AliEmcalJetTask.cxx
diff --git a/PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx b/PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx
new file mode 100644 (file)
index 0000000..86c8a74
--- /dev/null
@@ -0,0 +1,733 @@
+// $Id$
+//
+// Jet deltaPt task.
+//
+// Author: 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 "AliVParticle.h"
+#include "AliVTrack.h"
+#include "AliEmcalJet.h"
+#include "AliVEventHandler.h"
+#include "AliRhoParameter.h"
+#include "AliLog.h"
+
+#include "AliAnalysisTaskDeltaPt.h"
+
+ClassImp(AliAnalysisTaskDeltaPt)
+
+//________________________________________________________________________
+AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() : 
+  AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
+  fMCAna(kFALSE),
+  fMinRC2LJ(-1),
+  fEmbJetsName(""),
+  fEmbTracksName(""),
+  fEmbCaloName(""),
+  fRandTracksName("TracksRandomized"),
+  fRandCaloName("CaloClustersRandomized"),
+  fRCperEvent(-1),
+  fEmbJets(0),
+  fEmbTracks(0),
+  fEmbCaloClusters(0),
+  fRandTracks(0),
+  fRandCaloClusters(0),
+  fEmbeddedClusterId(-1),
+  fEmbeddedTrackId(-1),
+  fHistRCPhiEta(0),
+  fHistRCPtExLJVSDPhiLJ(0),
+  fHistRhoVSRCPt(0),
+  fHistEmbJetPhiEta(0),
+  fHistEmbPartPhiEta(0),
+  fHistRhoVSEmbBkg(0)
+{
+  // Default constructor.
+
+  for (Int_t i = 0; i < 4; i++) {
+    fHistRCPt[i] = 0;
+    fHistRCPtExLJ[i] = 0;
+    fHistRCPtRand[i] = 0; 
+    fHistDeltaPtRC[i] = 0;
+    fHistDeltaPtRCExLJ[i] = 0;
+    fHistDeltaPtRCRand[i] = 0;
+    fHistEmbNotFoundPhiEta[i] = 0;
+    fHistEmbJetsPt[i] = 0;
+    fHistEmbJetsCorrPt[i] = 0;
+    fHistEmbJetsArea[i] = 0;
+    fHistEmbPartPt[i] = 0;
+    fHistDistEmbPartJetAxis[i] = 0;
+    fHistDeltaPtEmb[i] = 0;
+  }
+
+  SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) : 
+  AliAnalysisTaskEmcalJet(name, kTRUE),
+  fMCAna(kFALSE),
+  fMinRC2LJ(-1),
+  fEmbJetsName(""),
+  fEmbTracksName(""),
+  fEmbCaloName(""),
+  fRandTracksName("TracksRandomized"),
+  fRandCaloName("CaloClustersRandomized"),
+  fRCperEvent(-1),
+  fEmbJets(0),
+  fEmbTracks(0),
+  fEmbCaloClusters(0),
+  fRandTracks(0),
+  fRandCaloClusters(0),
+  fEmbeddedClusterId(-1),
+  fEmbeddedTrackId(-1),
+  fHistRCPhiEta(0),
+  fHistRCPtExLJVSDPhiLJ(0),
+  fHistRhoVSRCPt(0),
+  fHistEmbJetPhiEta(0),
+  fHistEmbPartPhiEta(0),
+  fHistRhoVSEmbBkg(0)
+{
+  // Standard constructor.
+
+  for (Int_t i = 0; i < 4; i++) {
+    fHistRCPt[i] = 0;
+    fHistRCPtExLJ[i] = 0;
+    fHistRCPtRand[i] = 0; 
+    fHistDeltaPtRC[i] = 0;
+    fHistDeltaPtRCExLJ[i] = 0;
+    fHistDeltaPtRCRand[i] = 0;
+    fHistEmbNotFoundPhiEta[i] = 0;
+    fHistEmbJetsPt[i] = 0;
+    fHistEmbJetsCorrPt[i] = 0;
+    fHistEmbJetsArea[i] = 0;
+    fHistEmbPartPt[i] = 0;
+    fHistDistEmbPartJetAxis[i] = 0;
+    fHistDeltaPtEmb[i] = 0;
+  }
+
+  SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskDeltaPt::~AliAnalysisTaskDeltaPt()
+{
+  // Destructor.
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
+{
+  // Create user output.
+
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+  fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
+  fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
+  fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistRCPhiEta);
+
+  fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 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, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+  fHistRhoVSRCPt->GetXaxis()->SetTitle("A#rho [GeV/c]");
+  fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
+  fOutput->Add(fHistRhoVSRCPt);
+
+  if (!fJetsName.IsNull()) {
+    fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
+    fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
+    fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
+    fOutput->Add(fHistEmbJetPhiEta);
+    
+    fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
+    fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
+    fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
+    fOutput->Add(fHistEmbPartPhiEta);
+    
+    fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+    fHistRhoVSEmbBkg->GetXaxis()->SetTitle("A#rho [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 = "fHistRCPt_";
+    histname += i;
+    fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 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, fMinBinPt, fMaxBinPt * 2);
+    fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
+    fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistRCPtExLJ[i]);
+
+    histname = "fHistRCPtRand_";
+    histname += i;
+    fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
+    fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
+    fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistRCPtRand[i]);
+
+    histname = "fHistDeltaPtRC_";
+    histname += i;
+    fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
+    fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaPtRC[i]);
+
+    histname = "fHistDeltaPtRCExLJ_";
+    histname += i;
+    fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
+    fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaPtRCExLJ[i]);
+
+    histname = "fHistDeltaPtRCRand_";
+    histname += i;
+    fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
+    fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaPtRCRand[i]);
+
+    if (!fEmbJetsName.IsNull()) {
+      histname = "fHistEmbJetsPt_";
+      histname += i;
+      fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
+      fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
+      fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistEmbJetsPt[i]);
+
+      histname = "fHistEmbJetsCorrPt_";
+      histname += i;
+      fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+      fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
+      fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistEmbJetsCorrPt[i]);
+
+      histname = "fHistEmbJetsArea_";
+      histname += i;
+      fHistEmbJetsArea[i] = new TH1F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
+      fHistEmbJetsArea[i]->GetXaxis()->SetTitle("area");
+      fHistEmbJetsArea[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistEmbJetsArea[i]);
+
+      histname = "fHistEmbPartPt_";
+      histname += i;
+      fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
+      fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
+      fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistEmbPartPt[i]);
+
+      histname = "fHistDistEmbPartJetAxis_";
+      histname += i;
+      fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
+      fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
+      fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistDistEmbPartJetAxis[i]);
+
+      histname = "fHistEmbNotFoundPhiEta_";
+      histname += i;
+      fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
+      fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
+      fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
+      fOutput->Add(fHistEmbNotFoundPhiEta[i]);
+      
+      histname = "fHistDeltaPtEmb_";
+      histname += i;
+      fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+      fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [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
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
+{
+  // Fill histograms.
+
+  Int_t *sortedJets = GetSortedArray(fJets);
+  
+  AliEmcalJet* jet = 0;
+
+  if (sortedJets && sortedJets[0] > 0) 
+    jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
+
+  // ************
+  // Random cones
+  // _________________________________
+  
+  const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
+
+  // Simple random cones
+  for (Int_t i = 0; i < fRCperEvent; i++) {
+    Float_t RCpt = 0;
+    Float_t RCeta = 0;
+    Float_t RCphi = 0;
+    GetRandomCone(RCpt, RCeta, RCphi, 0);
+    if (RCpt > 0) {
+      fHistRCPhiEta->Fill(RCeta, RCphi);
+      fHistRhoVSRCPt->Fill(fRhoVal * rcArea, RCpt);
+
+      fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
+      fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
+    }
+  
+    // Random cones far from leading jet
+    RCpt = 0;
+    RCeta = 0;
+    RCphi = 0;
+    GetRandomCone(RCpt, RCeta, RCphi, jet);
+    if (RCpt > 0) {
+      if (jet) {
+       Float_t dphi = RCphi - jet->Phi();
+       if (dphi > 4.8) dphi -= TMath::Pi() * 2;
+       if (dphi < -1.6) dphi += TMath::Pi() * 2; 
+       fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
+      }
+      fHistRCPtExLJ[fCentBin]->Fill(RCpt / rcArea);
+      fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
+    }
+    
+    // Random cones with randomized particles
+    RCpt = 0;
+    RCeta = 0;
+    RCphi = 0;
+    GetRandomCone(RCpt, RCeta, RCphi, 0, fRandTracks, fRandCaloClusters);
+    if (RCpt > 0) {
+      fHistRCPtRand[fCentBin]->Fill(RCpt / rcArea);
+      fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
+    }  
+  }
+  
+  // ************
+  // Embedding
+  // _________________________________
+
+  if (!fJets)
+    return kTRUE;
+
+  AliEmcalJet *embJet  = 0;
+  TObject     *embPart = 0;
+
+  DoEmbJetLoop(embJet, embPart);
+
+  if (embJet) {
+    Double_t probePt = 0;
+    Double_t probeEta = 0;
+    Double_t probePhi = 0;
+
+    AliVCluster *cluster = dynamic_cast<AliVCluster*>(embPart);
+    if (cluster) {
+      TLorentzVector nPart;
+      cluster->GetMomentum(nPart, fVertex);
+
+      probeEta = nPart.Eta();
+      probePhi = nPart.Phi();
+      probePt = nPart.Pt();
+    }
+    else {
+      AliVTrack *track = dynamic_cast<AliVTrack*>(embPart);
+      if (track) {
+       probeEta = track->Eta();
+       probePhi = track->Phi();
+       probePt = track->Pt();
+      }
+      else {
+       AliWarning(Form("%s - Embedded jet found but embedded particle not found (wrong type?)!", GetName()));
+       return kTRUE;
+      }
+    }
+    Double_t distProbeJet = TMath::Sqrt((embJet->Eta() - probeEta) * (embJet->Eta() - probeEta) + (embJet->Phi() - probePhi) * (embJet->Phi() - probePhi));
+
+    fHistEmbPartPt[fCentBin]->Fill(probePt);
+    fHistEmbPartPhiEta->Fill(probeEta, probePhi);
+    
+    fHistDistEmbPartJetAxis[fCentBin]->Fill(distProbeJet);
+
+    fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
+    fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - fRhoVal * embJet->Area());
+    fHistEmbJetsArea[fCentBin]->Fill(embJet->Area());
+    fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
+
+    fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRhoVal - probePt);
+    fHistRhoVSEmbBkg->Fill(embJet->Area() * fRhoVal, embJet->Pt() - probePt);
+  }
+  else {
+    if (fEmbTracks)
+      DoEmbTrackLoop();
+    if (fEmbCaloClusters)
+      DoEmbClusterLoop();
+    if (fEmbeddedTrackId >= 0) {
+      AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
+      fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
+    }
+    else if (fEmbeddedClusterId >= 0) {
+      AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
+      TLorentzVector nPart;
+      cluster2->GetMomentum(nPart, fVertex);
+      fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
+    }
+    else {
+      AliWarning(Form("%s - Embedded particle not found!", GetName()));
+    }
+  }
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
+{
+  // Do track loop.
+
+  if (!fEmbTracks)
+    return;
+
+  fEmbeddedTrackId = -1;
+
+  Int_t ntracks = fEmbTracks->GetEntriesFast();
+
+  for (Int_t i = 0; i < ntracks; i++) {
+
+    AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track  
+
+    if (!track) {
+      AliError(Form("Could not retrieve track %d",i)); 
+      continue; 
+    }
+
+    AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
+    
+    if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
+      continue;
+
+    if (track->GetLabel() == 100) {
+      fEmbeddedTrackId = i;
+      continue;
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
+{
+  // Do cluster loop.
+
+  if (!fEmbCaloClusters)
+    return;
+
+  fEmbeddedClusterId = -1;
+
+  Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
+
+  for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
+    AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
+    if (!cluster) {
+      AliError(Form("Could not receive cluster %d", iClusters));
+      continue;
+    }  
+
+    if (!AcceptCluster(cluster, kTRUE)) 
+      continue;
+
+    if (cluster->Chi2() == 100) {
+      fEmbeddedClusterId = iClusters;
+      continue;
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
+{
+  // Do the embedded jet loop.
+
+  if (!fEmbJets)
+    return;
+
+  TLorentzVector maxClusVect;
+
+  Int_t nembjets = fEmbJets->GetEntriesFast();
+
+  for (Int_t ij = 0; ij < nembjets; ij++) {
+      
+    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
+      
+    if (!jet) {
+      AliError(Form("Could not receive jet %d", ij));
+      continue;
+    } 
+      
+    if (!AcceptJet(jet))
+      continue;
+
+    if (!jet->IsMC())
+      continue;
+
+    AliVParticle *maxTrack = 0;
+
+    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
+      AliVParticle *track = jet->TrackAt(it, fEmbTracks);
+      
+      if (!track) 
+       continue;
+
+      if (track->GetLabel() != 100)
+       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, fEmbCaloClusters);
+      
+      if (!cluster) 
+       continue;
+
+      if (cluster->Chi2() != 100)
+       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)) {
+      embPart = maxClus;
+      embJet = jet;
+    }
+    else if (maxTrack) {
+      embPart = maxTrack;
+      embJet = jet;
+    }
+
+    return;  // MC jet found, exit
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
+                                       AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
+{
+  // Get rigid cone.
+
+  if (!tracks)
+    tracks = fTracks;
+
+  if (!clusters)
+    clusters = fCaloClusters;
+
+  eta = 0;
+  phi = 0;
+  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 = fJetMaxEta;
+  Float_t minEta = fJetMinEta;
+  Float_t maxPhi = fJetMaxPhi;
+  Float_t minPhi = fJetMinPhi;
+
+  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) {
+    AliWarning(Form("%s: Could not get random cone!", GetName()));
+    return;
+  }
+
+  if (clusters) {
+    Int_t nclusters =  clusters->GetEntriesFast();
+    for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
+      AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
+      if (!cluster) {
+       AliError(Form("Could not receive cluster %d", iClusters));
+       continue;
+      }  
+      
+      if (!AcceptCluster(cluster, fMCAna))
+       continue;
+      
+      TLorentzVector nPart;
+      cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
+     
+      Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
+
+      if (d <= fJetRadius) 
+       pt += nPart.Pt();
+    }
+  }
+
+  if (tracks) {
+    Int_t ntracks = tracks->GetEntriesFast();
+    for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
+      AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));         
+      if(!track) {
+       AliError(Form("Could not retrieve track %d",iTracks)); 
+       continue; 
+      }
+      
+      if (!AcceptTrack(track, fMCAna)) 
+       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();
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::ExecOnce()
+{
+  // Initialize the analysis.
+
+  if (!fEmbJetsName.IsNull() && !fEmbJets) {
+    fEmbJets =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
+    if (!fEmbJets) {
+      AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
+      return;
+    }
+    else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
+      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data())); 
+      fEmbJets = 0;
+      return;
+    }
+  }
+
+  if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
+    fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
+    if (!fEmbCaloClusters) {
+      AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
+      return;
+    }
+    else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
+      AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data())); 
+      fEmbCaloClusters = 0;
+      return;
+    }
+  }
+
+  if (!fEmbTracksName.IsNull() && !fEmbTracks) {
+    fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
+    if (!fEmbTracks) {
+      AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
+      return;
+    }
+    else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
+      AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data())); 
+      fEmbTracks = 0;
+      return;
+    }
+  }
+
+  if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
+    fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
+    if (!fRandCaloClusters) {
+      AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
+      return;
+    }
+    else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
+      AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data())); 
+      fRandCaloClusters = 0;
+      return;
+    }
+  }
+
+  if (!fRandTracksName.IsNull() && !fRandTracks) {
+    fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
+    if (!fRandTracks) {
+      AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
+      return;
+    }
+    else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
+      AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data())); 
+      fRandTracks = 0;
+      return;
+    }
+  }
+
+  AliAnalysisTaskEmcalJet::ExecOnce();
+
+  if (fRCperEvent < 0) {
+    Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
+    Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
+    fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
+    if (fRCperEvent == 0)
+      fRCperEvent = 1;
+  }
+
+  if (fMinRC2LJ < 0)
+    fMinRC2LJ = fJetRadius * 1.5;
+
+  const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 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;
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::Terminate(Option_t *) 
+{
+  // Called once at the end of the analysis.
+}
diff --git a/PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.h b/PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.h
new file mode 100644 (file)
index 0000000..b344d01
--- /dev/null
@@ -0,0 +1,89 @@
+#ifndef ALIANALYSISTASKDELTAPT_H
+#define ALIANALYSISTASKDELTAPT_H
+
+// $Id$
+
+class TClonesArray;
+class TString;
+class TH1F;
+class TH2F;
+class AliRhoParameter;
+
+#include "AliAnalysisTaskEmcalJet.h"
+
+class AliAnalysisTaskDeltaPt : public AliAnalysisTaskEmcalJet {
+ public:
+
+  AliAnalysisTaskDeltaPt();
+  AliAnalysisTaskDeltaPt(const char *name);
+  virtual ~AliAnalysisTaskDeltaPt();
+
+  void                        UserCreateOutputObjects();
+  void                        Terminate(Option_t *option);
+
+  void                        SetJetMinRC2LJ(Float_t d)                            { fMinRC2LJ                = d          ; }
+  void                        SetEmbJetsName(const char *n)                        { fEmbJetsName             = n          ; } 
+  void                        SetEmbTracksName(const char *n)                      { fEmbTracksName           = n          ; }
+  void                        SetEmbClusName(const char *n)                        { fEmbCaloName             = n          ; }
+  void                        SetRandTracksName(const char *n)                     { fRandTracksName          = n          ; }
+  void                        SetRandClusName(const char *n)                       { fRandCaloName            = n          ; }
+  void                        SetMC(Bool_t m)                                      { fMCAna                   = m          ; }
+  void                        SetRCperEvent(Int_t n)                               { fRCperEvent              = n          ; }
+
+ protected:
+  void                        ExecOnce()                                                                                    ;
+  Bool_t                      FillHistograms()                                                                              ;
+  void                        GetLeadingJets(Int_t &maxJetIndex, Int_t &max2JetIndex)                                       ;
+  void                        DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)                                            ;
+  void                        DoEmbTrackLoop()                                                                              ;
+  void                        DoEmbClusterLoop()                                                                            ;
+  void                        GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi, 
+                                           AliEmcalJet *jet = 0, TClonesArray* tracks = 0, TClonesArray* clusters = 0) const;
+
+  Bool_t                      fMCAna;                      // =true MC analysis (toy model)
+  Float_t                     fMinRC2LJ;                   // Minimum distance random cone to leading jet
+  TString                     fEmbJetsName;                // Name of embedded jet collection
+  TString                     fEmbTracksName;              // Name of embedded track collection
+  TString                     fEmbCaloName;                // Name of embedded calo cluster collection
+  TString                     fRandTracksName;             // Name of randomized track collection
+  TString                     fRandCaloName;               // Name of randomized calo cluster collection
+  Int_t                       fRCperEvent;                 // No. of random cones per event
+
+  TClonesArray               *fEmbJets;                    //!Embedded jets
+  TClonesArray               *fEmbTracks;                  //!Embedded tracks
+  TClonesArray               *fEmbCaloClusters;            //!Embedded clusters  
+  TClonesArray               *fRandTracks;                 //!Randomized tracks
+  TClonesArray               *fRandCaloClusters;           //!Randomized clusters
+  Int_t                       fEmbeddedClusterId;          //!Embedded cluster id
+  Int_t                       fEmbeddedTrackId;            //!Embedded track id
+
+  // Random cones
+  TH2F                       *fHistRCPhiEta;               //!Phi-Eta distribution of random cones
+  TH1F                       *fHistRCPt[4];                //!Random cone pt
+  TH1F                       *fHistRCPtExLJ[4];            //!Random cone pt, imposing min distance from leading jet
+  TH1F                       *fHistRCPtRand[4];            //!Random cone pt, randomized particles
+  TH2F                       *fHistRCPtExLJVSDPhiLJ;       //!Random cone pt, imposing min distance from leading jet, vs. deltaPhi leading jet
+  TH2F                       *fHistRhoVSRCPt;              //!Area(RC) * rho vs. Pt(RC)
+  TH1F                       *fHistDeltaPtRC[4];           //!deltaPt = Pt(RC) - A * rho
+  TH1F                       *fHistDeltaPtRCExLJ[4];       //!deltaPt = Pt(RC) - A * rho, imposing min distance from leading jet
+  TH1F                       *fHistDeltaPtRCRand[4];       //!deltaPt = Pt(RC) - A * rho, randomzied particles
+
+  // Jet embedding
+  TH2F                       *fHistEmbNotFoundPhiEta[4];   //!Phi-Eta of "not found" embedded particles
+  TH1F                       *fHistEmbJetsPt[4];           //!Pt distribution of embedded jets
+  TH1F                       *fHistEmbJetsCorrPt[4];       //!Pt-rho*A distribution of embedded jets
+  TH1F                       *fHistEmbJetsArea[4];         //!Area distribution of embedded jets
+  TH1F                       *fHistEmbPartPt[4];           //!Pt distribution of embedded particle
+  TH2F                       *fHistEmbJetPhiEta;           //!Phi-Eta distribution of embedded jets
+  TH2F                       *fHistEmbPartPhiEta;          //!Phi-Eta distribution of embedded particles
+  TH1F                       *fHistDistEmbPartJetAxis[4];  //!Distance between embedded particle and jet axis
+  TH2F                       *fHistRhoVSEmbBkg;            //!Area(embjet) * rho vs. Pt(embjet) - Pt(embtrack)
+  TH1F                       *fHistDeltaPtEmb[4];          //!deltaPt = Pt(embjet) - Area(embjet) * rho - Pt(embtrack)
+
+ private:
+  AliAnalysisTaskDeltaPt(const AliAnalysisTaskDeltaPt&);            // not implemented
+  AliAnalysisTaskDeltaPt &operator=(const AliAnalysisTaskDeltaPt&); // not implemented
+
+  ClassDef(AliAnalysisTaskDeltaPt, 1) // deltaPt analysis task
+};
+#endif
index 364a6215bcac699f10db113c875c59c3bbe4f40a..a9a831101889e2bfc42d8c59a2b3c66d2125f8a8 100644 (file)
@@ -29,18 +29,22 @@ AliAnalysisTaskEmcalJet::AliAnalysisTaskEmcalJet() :
   AliAnalysisTaskEmcal("AliAnalysisTaskEmcalJet"),
   fJetRadius(0.4),
   fJetsName(),
-  fPtBiasJetTrack(5),
-  fPtBiasJetClus(5),
+  fRhoName(),
+  fPtBiasJetTrack(0),
+  fPtBiasJetClus(0),
   fJetPtCut(1),
-  fJetAreaCut(-1),
-  fPercAreaCut(0.8),
-  fMinEta(-0.9),
-  fMaxEta(0.9),
-  fMinPhi(-10),
-  fMaxPhi(10),
+  fJetAreaCut(0.4),
+  fPercAreaCut(-1),
+  fAreaEmcCut(0),
+  fJetMinEta(-0.9),
+  fJetMaxEta(0.9),
+  fJetMinPhi(-10),
+  fJetMaxPhi(10),
   fMaxClusterPt(100),
   fMaxTrackPt(100),
-  fJets(0)
+  fJets(0),
+  fRho(0),
+  fRhoVal(0)
 {
   // Default constructor.
 }
@@ -50,18 +54,22 @@ AliAnalysisTaskEmcalJet::AliAnalysisTaskEmcalJet(const char *name, Bool_t histo)
   AliAnalysisTaskEmcal(name, histo),
   fJetRadius(0.4),
   fJetsName(),
-  fPtBiasJetTrack(5),
-  fPtBiasJetClus(5),
+  fRhoName(),
+  fPtBiasJetTrack(0),
+  fPtBiasJetClus(0),
   fJetPtCut(1),
-  fJetAreaCut(-1),
-  fPercAreaCut(0.8),
-  fMinEta(-0.9),
-  fMaxEta(0.9),
-  fMinPhi(-10),
-  fMaxPhi(10),
+  fJetAreaCut(0.4),
+  fPercAreaCut(-1),
+  fAreaEmcCut(0),
+  fJetMinEta(-0.9),
+  fJetMaxEta(0.9),
+  fJetMinPhi(-10),
+  fJetMaxPhi(10),
   fMaxClusterPt(100),
   fMaxTrackPt(100),
-  fJets(0)
+  fJets(0),
+  fRho(0),
+  fRhoVal(0)
 {
   // Standard constructor.
 }
@@ -77,14 +85,14 @@ Bool_t AliAnalysisTaskEmcalJet::AcceptBiasJet(AliEmcalJet *jet) const
 { 
   // Accept jet with a bias.
 
-  if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet->MaxClusterPt() < fPtBiasJetClus))
+  if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus)
     return kFALSE;
   else
     return kTRUE;
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskEmcalJet::AcceptJet(AliEmcalJet *jet, Bool_t bias, Bool_t upCut) const
+Bool_t AliAnalysisTaskEmcalJet::AcceptJet(AliEmcalJet *jet) const
 {   
   // Return true if jet is accepted.
 
@@ -92,12 +100,14 @@ Bool_t AliAnalysisTaskEmcalJet::AcceptJet(AliEmcalJet *jet, Bool_t bias, Bool_t
     return kFALSE;
   if (jet->Area() <= fJetAreaCut)
     return kFALSE;
-  if (bias && !AcceptBiasJet(jet))
+  if (jet->AreaEmc()<fAreaEmcCut)
     return kFALSE;
-  if (upCut && (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt))
+  if (!AcceptBiasJet(jet))
+    return kFALSE;
+  if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
     return kFALSE;
 
-  return (Bool_t)(jet->Eta() > fMinEta && jet->Eta() < fMaxEta && jet->Phi() > fMinPhi && jet->Phi() < fMaxPhi);
+  return (Bool_t)(jet->Eta() > fJetMinEta && jet->Eta() < fJetMaxEta && jet->Phi() > fJetMinPhi && jet->Phi() < fJetMaxPhi);
 }
 
 //________________________________________________________________________
@@ -122,6 +132,8 @@ void AliAnalysisTaskEmcalJet::ExecOnce()
 {
   // Init the analysis.
 
+  AliAnalysisTaskEmcal::ExecOnce();
+
   if (fPercAreaCut >= 0) {
     if (fJetAreaCut >= 0)
       AliInfo(Form("%s: jet area cut will be calculated as a percentage of the average area, given value will be overwritten", GetName()));
@@ -129,28 +141,125 @@ void AliAnalysisTaskEmcalJet::ExecOnce()
   }
 
   if (fAnaType == kTPC) {
-    SetEtaLimits(-0.9 + fJetRadius, 0.9 - fJetRadius);
-    SetPhiLimits(-10, 10);
-  } else if (fAnaType == kEMCAL || fAnaType == kTPCSmall || fAnaType == kEMCALOnly) {
-    AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
-    if (geom) {
-      SetEtaLimits(geom->GetArm1EtaMin() + fJetRadius, geom->GetArm1EtaMax() - fJetRadius);
-      SetPhiLimits(geom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, geom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
+    SetJetEtaLimits(-0.5, 0.5);
+    SetJetPhiLimits(-10, 10);
+    SetTrackEtaLimits(-0.5 - fJetRadius, 0.5 + fJetRadius);
+    SetTrackPhiLimits(-10, 10);
+  } 
+  else if (fAnaType == kEMCAL && fGeom) {
+    SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
+    SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
+    SetTrackEtaLimits(fGeom->GetArm1EtaMin(), fGeom->GetArm1EtaMax());
+    SetTrackPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad(), fGeom->GetArm1PhiMax() * TMath::DegToRad());
+  }
+
+  if (!fRhoName.IsNull() && !fRho) {
+    fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
+    if (!fRho) {
+      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
+      fInitialized = kFALSE;
+      return;
+    }
+  }
+
+  if (!fJetsName.IsNull() && !fJets) {
+    fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
+    if (!fJets) {
+      AliError(Form("%s: Could not retrieve jets %s!", GetName(), fJetsName.Data()));
+      fInitialized = kFALSE;
+      return;
     }
-    else {
-      AliWarning(Form("%s: Can not create geometry", GetName()));
+    else if (!fJets->GetClass()->GetBaseClass("AliEmcalJet")) {
+      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetsName.Data())); 
+      fJets = 0;
+      fInitialized = kFALSE;
+      return;
     }
-  } else {
-    AliWarning(Form("%s: Analysis type not recognized! Assuming kTPC!", GetName()));
-    SetAnaType(kTPC);
-    ExecOnce();
-    return;
   }
+}
 
-  if (fAnaType == kTPCSmall)
-    fAnaType = kTPC;
+//________________________________________________________________________
+Int_t* AliAnalysisTaskEmcalJet::GetSortedArray(TClonesArray *array) const
+{
+  // Get the leading jets.
 
-  AliAnalysisTaskEmcal::ExecOnce();
+  static Float_t pt[9999];
+  static Int_t indexes[9999];
+
+  if (!array)
+    return 0;
+
+  const Int_t n = array->GetEntriesFast();
+
+  if (fJets->GetClass()->GetBaseClass("AliEmcalJet")) {
+
+    for (Int_t i = 0; i < n; i++) {
+
+      pt[i] = -FLT_MAX;
+
+      AliEmcalJet* jet = static_cast<AliEmcalJet*>(array->At(i));
+      
+      if (!jet) {
+       AliError(Form("Could not receive jet %d", i));
+       continue;
+      }
+      
+      if (!AcceptJet(jet))
+       continue;
+      
+      pt[i] = jet->Pt() - fRhoVal * jet->Area();
+    }
+  }
+
+  else if (fJets->GetClass()->GetBaseClass("AliVTrack")) {
+
+    for (Int_t i = 0; i < n; i++) {
+
+      pt[i] = -FLT_MAX;
+
+      AliVTrack* track = static_cast<AliVTrack*>(array->At(i));
+      
+      if (!track) {
+       AliError(Form("Could not receive track %d", i));
+       continue;
+      }  
+      
+      if (!AcceptTrack(track))
+       continue;
+      
+      pt[i] = track->Pt();
+    }
+  }
+
+  else if (fJets->GetClass()->GetBaseClass("AliVCluster")) {
+
+    for (Int_t i = 0; i < n; i++) {
+
+      pt[i] = -FLT_MAX;
+
+      AliVCluster* cluster = static_cast<AliVCluster*>(array->At(i));
+      
+      if (!cluster) {
+       AliError(Form("Could not receive cluster %d", i));
+       continue;
+      }  
+      
+      if (!AcceptCluster(cluster))
+       continue;
+
+      TLorentzVector nPart;
+      cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
+      
+      pt[i] = nPart.Pt();
+    }
+  }
+
+  TMath::Sort(n, pt, indexes);
+
+  if (pt[indexes[0]] == -FLT_MAX) 
+    return 0;
+
+  return indexes;
 }
 
 //________________________________________________________________________
@@ -191,18 +300,8 @@ Bool_t AliAnalysisTaskEmcalJet::RetrieveEventObjects()
   if (!AliAnalysisTaskEmcal::RetrieveEventObjects())
     return kFALSE;
 
-  if (!fJetsName.IsNull() && !fJets) {
-    fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
-    if (!fJets) {
-      AliError(Form("%s: Could not retrieve jets %s!", GetName(), fJetsName.Data()));
-      return kFALSE;
-    }
-    else if (!fJets->GetClass()->GetBaseClass("AliEmcalJet")) {
-      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetsName.Data())); 
-      fJets = 0;
-      return kFALSE;
-    }
-  }
+  if (fRho)
+    fRhoVal = fRho->GetVal();
 
   return kTRUE;
 }
index 26cd04f86323b8324b71cec706f81bfe27a05cda..f817e4cb7fadf195d1422aaf813d961c6b83235d 100644 (file)
@@ -19,47 +19,53 @@ class AliAnalysisTaskEmcalJet : public AliAnalysisTaskEmcal {
   AliAnalysisTaskEmcalJet(const char *name, Bool_t histo=kFALSE); 
   virtual ~AliAnalysisTaskEmcalJet();
 
-  void                        SetEtaLimits(Float_t min, Float_t max)               { fMinEta = min, fMaxEta = max ; }
-  void                        SetJetAreaCut(Float_t cut)                           { fJetAreaCut     = cut        ; }
-  void                        SetPercAreaCut(Float_t p)                            { fPercAreaCut    = p          ; }
-  void                        SetJetPtCut(Float_t cut)                             { fJetPtCut       = cut        ; }
-  void                        SetJetRadius(Float_t r)                              { fJetRadius      = r          ; } 
-  void                        SetJetsName(const char *n)                           { fJetsName       = n          ; }
-  void                        SetMaxClusterPt(Float_t b)                           { fMaxClusterPt  = b           ; }
-  void                        SetMaxTrackPt(Float_t b)                             { fMaxTrackPt = b              ; }
-  void                        SetPhiLimits(Float_t min, Float_t max)               { fMinPhi = min, fMaxPhi = max ; }
-  void                        SetPtBiasJetClus(Float_t b)                          { fPtBiasJetClus  = b          ; }
-  void                        SetPtBiasJetTrack(Float_t b)                         { fPtBiasJetTrack = b          ; }
-
+  void                        SetJetEtaLimits(Float_t min, Float_t max)            { fJetMinEta = min, fJetMaxEta = max ; }
+  void                        SetJetPhiLimits(Float_t min, Float_t max)            { fJetMinPhi = min, fJetMaxPhi = max ; }
+  void                        SetJetAreaCut(Float_t cut)                           { fJetAreaCut     = cut              ; }
+  void                        SetPercAreaCut(Float_t p)                            { fPercAreaCut    = p                ; }
+  void                        SetAreaEmcCut(Double_t a = 0.99)                     { fAreaEmcCut     = a                ; }
+  void                        SetJetPtCut(Float_t cut)                             { fJetPtCut       = cut              ; }
+  void                        SetJetRadius(Float_t r)                              { fJetRadius      = r                ; } 
+  void                        SetJetsName(const char *n)                           { fJetsName       = n                ; }
+  virtual void                SetRhoName(const char *n)                            { fRhoName        = n                ; }
+  void                        SetMaxClusterPt(Float_t b)                           { fMaxClusterPt   = b                ; }
+  void                        SetMaxTrackPt(Float_t b)                             { fMaxTrackPt     = b                ; }
+  void                        SetPtBiasJetClus(Float_t b)                          { fPtBiasJetClus  = b                ; }
+  void                        SetPtBiasJetTrack(Float_t b)                         { fPtBiasJetTrack = b                ; }
  
  protected:
-  virtual Bool_t              AcceptJet(AliEmcalJet* jet, Bool_t bias = kTRUE, Bool_t upCut = kTRUE)   const;
+  virtual Bool_t              AcceptJet(AliEmcalJet* jet)                                              const;
   Bool_t                      AcceptBiasJet(AliEmcalJet* jet)                                          const;
   void                        ExecOnce()                                                                    ;
-  AliRhoParameter            *GetRhoFromEvent(const char *name);
+  AliRhoParameter            *GetRhoFromEvent(const char *name)                                             ;
+  Int_t                      *GetSortedArray(TClonesArray *array)                                      const;
   Bool_t                      IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted = kTRUE)        const;
   Bool_t                      IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted = kTRUE)       const;
   Bool_t                      RetrieveEventObjects()                                                        ;
 
   Float_t                     fJetRadius;                  // jet radius
   TString                     fJetsName;                   // name of jet collection
+  TString                     fRhoName;                    // Name of rho object
   Float_t                     fPtBiasJetTrack;             // select jets with a minimum pt track
   Float_t                     fPtBiasJetClus;              // select jets with a minimum pt cluster
   Float_t                     fJetPtCut;                   // cut on jet pt
   Float_t                     fJetAreaCut;                 // cut on jet area
   Float_t                     fPercAreaCut;                // cut on jet area as a percentage of average jet area
-  Float_t                     fMinEta;                     // minimum eta jet acceptance
-  Float_t                     fMaxEta;                     // maximum eta jet acceptance
-  Float_t                     fMinPhi;                     // minimum phi jet acceptance
-  Float_t                     fMaxPhi;                     // maximum phi jet acceptance  
+  Float_t                     fAreaEmcCut;                 // minimum cut on jet emcal area
+  Float_t                     fJetMinEta;                  // minimum eta jet acceptance
+  Float_t                     fJetMaxEta;                  // maximum eta jet acceptance
+  Float_t                     fJetMinPhi;                  // minimum phi jet acceptance
+  Float_t                     fJetMaxPhi;                  // maximum phi jet acceptance  
   Float_t                     fMaxClusterPt;               // maximum cluster constituent pt to accept the jet
   Float_t                     fMaxTrackPt;                 // maximum track constituent pt to accept the jet
   TClonesArray               *fJets;                       //!jets
+  AliRhoParameter            *fRho;                        //!Event rho
+  Double_t                    fRhoVal;                     //!Event rho value
 
  private:
   AliAnalysisTaskEmcalJet(const AliAnalysisTaskEmcalJet&);            // not implemented
   AliAnalysisTaskEmcalJet &operator=(const AliAnalysisTaskEmcalJet&); // not implemented
 
-  ClassDef(AliAnalysisTaskEmcalJet, 3) // EMCAL Jet base analysis task
+  ClassDef(AliAnalysisTaskEmcalJet, 4) // EMCAL Jet base analysis task
 };
 #endif
index 34e1caf6f52de0054b11026852fa052dcada803d..45263c4d3908132742372832445b06a3a342f848 100644 (file)
 #include "AliAnalysisTaskRho.h"
 
 #include <TClonesArray.h>
-#include <TF1.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TList.h>
 #include <TMath.h>
 
 #include "AliAnalysisManager.h"
-#include "AliCentrality.h"
 #include "AliEmcalJet.h"
 #include "AliLog.h"
 #include "AliRhoParameter.h"
-#include "AliVCluster.h"
-#include "AliVEventHandler.h"
 
 ClassImp(AliAnalysisTaskRho)
 
 //________________________________________________________________________
 AliAnalysisTaskRho::AliAnalysisTaskRho() : 
-  AliAnalysisTaskRhoBase(),
-  fTracksName(),
-  fJetsName(),
-  fRhoScaledName(""),
-  fPhiMin(0),
-  fPhiMax(0),
-  fEtaMin(0),
-  fEtaMax(0),
-  fAreaCut(0),
-  fAreaEmcCut(0),
-  fNExclLeadJets(0),
-  fScaleFunction(0),
-  fCreateHisto(kFALSE),
-  fTracks(0),
-  fJets(0),
-  fOutputList(0),
-  fHistCentrality(0),
-  fHistJetPt(0),
-  fHistJetArea(0),
-  fHistRhovsCent(0),
-  fHistDeltaRhovsCent(0),
-  fHistDeltaRhoScalevsCent(0),
-  fHistJetPtvsCent(0),
-  fHistJetAreavsCent(0),
-  fHistNjetvsCent(0),
-  fHistRhovsNtrack(0),
-  fHistDeltaRhovsNtrack(0),
-  fHistDeltaRhoScalevsNtrack(0),
-  fHistJetPtvsNtrack(0),
-  fHistJetAreavsNtrack(0),
-  fHistNjetvsNtrack(0),
-  fRhoScaled(0)
+  AliAnalysisTaskRhoBase("AliAnalysisTaskRho"),
+  fNExclLeadJets(0)
 {
   // Constructor.
 }
 
 //________________________________________________________________________
 AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name, Bool_t histo) :
-  AliAnalysisTaskRhoBase(name),
-  fTracksName("tracks"),
-  fJetsName("KtJets"),
-  fRhoScaledName(""),
-  fPhiMin(0),
-  fPhiMax(TMath::TwoPi()),
-  fEtaMin(-0.5),
-  fEtaMax(+0.5),
-  fAreaCut(0.01),
-  fAreaEmcCut(0),
-  fNExclLeadJets(1),
-  fScaleFunction(0),
-  fCreateHisto(histo),
-  fTracks(0),
-  fJets(0),
-  fOutputList(0),
-  fHistCentrality(0),
-  fHistJetPt(0),
-  fHistJetArea(0),
-  fHistRhovsCent(0),
-  fHistDeltaRhovsCent(0),
-  fHistDeltaRhoScalevsCent(0),
-  fHistJetPtvsCent(0),
-  fHistJetAreavsCent(0),
-  fHistNjetvsCent(0),
-  fHistRhovsNtrack(0),
-  fHistDeltaRhovsNtrack(0),
-  fHistDeltaRhoScalevsNtrack(0),
-  fHistJetPtvsNtrack(0),
-  fHistJetAreavsNtrack(0),
-  fHistNjetvsNtrack(0),
-  fRhoScaled(0)
+  AliAnalysisTaskRhoBase(name, histo),
+  fNExclLeadJets(0)
 {
   // Constructor.
-
-  if (fCreateHisto)
-    DefineOutput(1, TList::Class());
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskRho::UserCreateOutputObjects()
+Bool_t AliAnalysisTaskRho::Run() 
 {
-  // User create output objects, called at the beginning of the analysis.
-
-  AliAnalysisTaskRhoBase::UserCreateOutputObjects();
-
-  if (fScaleFunction) {
-    fRhoScaledName = fRhoName;
-    fRhoScaledName += "_Scaled";
-    fRhoScaled = new AliRhoParameter(fRhoScaledName, 0);  
-  }
-
-  if (!fCreateHisto)
-    return;
-
-  OpenFile(1);
-  fOutputList = new TList();
-  fOutputList->SetOwner();
-
-  fHistCentrality             = new TH1F("Centrality",            "Centrality",            101, -1,  100);
-  fHistRhovsCent              = new TH2F("RhovsCent",             "RhovsCent",             101, -1,  100,   500,  0,   500);
-  fHistDeltaRhovsCent         = new TH2F("DeltaRhovsCent",        "DetlaRhovsCent",        101, -1,  100,   500, -250, 250);
-  fHistDeltaRhoScalevsCent    = new TH2F("DeltaRhoScalevsCent",   "DeltaRhoScalevsCent",   101, -1,  100,   500, -250, 250);
-  fHistJetPtvsCent            = new TH2F("JetPtvsCent",           "JetPtvsCent",           101, -1,  100,   200,  0,   500);
-  fHistJetAreavsCent          = new TH2F("JetAreavsCent",         "JetAreavsCent",         101, -1,  100,   100,  0,   1.0);
-  fHistNjetvsCent             = new TH2F("NjetvsCent",            "NjetvsCent",            101, -1,  100,   100,  0,   100);
-
-  fHistRhovsNtrack            = new TH2F("RhovsNtrack",           "RhovsNtrack",           500,  0,  2500,  500,  0,   500);
-  fHistDeltaRhovsNtrack       = new TH2F("DeltaRhovsNtrack",      "DeltaRhovsNtrack",      500,  0,  2500,  500, -250, 250);
-  fHistDeltaRhoScalevsNtrack  = new TH2F("DeltaRhoScalevsNtrack", "DeltaRhoScalevsNtrack", 500,  0,  2500,  500, -250, 250);
-  fHistJetPtvsNtrack          = new TH2F("JetPtvsNtrack",         "JetPtvsNtrack",         500,  0,  2500,  200,  0,   500);
-  fHistJetAreavsNtrack        = new TH2F("JetAreavsNtrack",       "JetAreavsNtrack",       500,  0,  2500,  100,  0,   1.0);
-  fHistNjetvsNtrack           = new TH2F("NjetvsNtrack",          "rNjetvsNtrack",         500,  0,  2500,  100,  0,   100);
-
-  fHistJetPt                  = new TH1F("JetPt",                  "Jet Pt",               100,   0,    250);
-  fHistJetArea                = new TH1F("JetArea",                "Jet Area",             100, 0.0,    1.0);
-  
-  fOutputList->Add(fHistCentrality);
-  fOutputList->Add(fHistRhovsCent);
-  fOutputList->Add(fHistJetPt);
-  fOutputList->Add(fHistJetArea);
-  fOutputList->Add(fHistDeltaRhovsCent);
-  fOutputList->Add(fHistDeltaRhoScalevsCent);
-  fOutputList->Add(fHistJetPtvsCent);
-  fOutputList->Add(fHistJetAreavsCent);
-  fOutputList->Add(fHistNjetvsCent);
-  
-  fOutputList->Add(fHistRhovsNtrack);
-  fOutputList->Add(fHistDeltaRhovsNtrack);
-  fOutputList->Add(fHistDeltaRhoScalevsNtrack);
-  fOutputList->Add(fHistJetPtvsNtrack);
-  fOutputList->Add(fHistJetAreavsNtrack);
-  fOutputList->Add(fHistNjetvsNtrack);
-  
-  PostData(1, fOutputList);
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskRho::UserExec(Option_t *) 
-{
-  // Main loop, called for each event.
-
-  if (!fIsInit) {
-    ExecOnce();
-    fIsInit = 1;
-  }
+  // Run the analysis.
 
   fRho->SetVal(0);
   if (fRhoScaled)
     fRhoScaled->SetVal(0);
 
   if (!fJets)
-    return;
-
-  DetermineCent();
+    return kFALSE;
 
   const Int_t Njets   = fJets->GetEntries();
 
@@ -194,13 +59,7 @@ void AliAnalysisTaskRho::UserExec(Option_t *)
        continue;
       } 
 
-      if (jet->Area() < fAreaCut)
-        continue;
-      if (jet->AreaEmc()<fAreaEmcCut)
-        continue;
-      if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
-        continue;
-      if ((jet->Eta() < fEtaMin) || (jet->Eta() > fEtaMax))
+      if (!AcceptJet(jet))
         continue;
 
       if (jet->Pt() > maxJetPts[0]) {
@@ -215,7 +74,7 @@ void AliAnalysisTaskRho::UserExec(Option_t *)
     }
     if (fNExclLeadJets < 2) {
       maxJetIds[1] = -1;
-      maxJetPts[1] = -1;
+      maxJetPts[1] = 0;
     }
   }
 
@@ -235,101 +94,23 @@ void AliAnalysisTaskRho::UserExec(Option_t *)
       continue;
     } 
 
-    if (jet->Area() < fAreaCut)
-      continue;
-    if (jet->AreaEmc()<fAreaEmcCut)
-      continue;
-    if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
-      continue;
-    if ((jet->Eta() < fEtaMin) || (jet->Eta() > fEtaMax))
+    if (!AcceptJet(jet))
       continue;
 
     rhovec[NjetAcc] = jet->Pt() / jet->Area();
     ++NjetAcc;
-
-    if (fCreateHisto) {
-      // filling histograms
-      const Int_t Ntracks = fTracks->GetEntriesFast();
-      fHistJetPt->Fill(jet->Pt());
-      fHistJetArea->Fill(jet->Area());
-      fHistJetPtvsCent->Fill(fCent, jet->Pt());
-      fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
-      fHistJetAreavsCent->Fill(fCent, jet->Area());
-      fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
-    }
   }
 
   if (NjetAcc > 0) {
     //find median value
-    Double_t rho0      = TMath::Median(NjetAcc, rhovec);
-    fRho->SetVal(rho0);
-    Double_t rhoScaled = rho0;
+    Double_t rho = TMath::Median(NjetAcc, rhovec);
+    fRho->SetVal(rho);
+
     if (fRhoScaled) {
-      rhoScaled *= GetScaleFactor(fCent);
+      Double_t rhoScaled = rho * GetScaleFactor(fCent);
       fRhoScaled->SetVal(rhoScaled);
     }
-
-    if (fCreateHisto) {
-      // filling other histograms
-      Double_t rho        = GetRhoFactor(fCent);
-      const Int_t Ntracks = fTracks->GetEntriesFast();
-      fHistRhovsCent->Fill(fCent, rho0);
-      fHistDeltaRhovsCent->Fill(fCent, rho0 - rho);
-      fHistDeltaRhoScalevsCent->Fill(fCent, rhoScaled - rho);
-      fHistRhovsNtrack->Fill(Ntracks, rho0);
-      fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rho);
-      fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rhoScaled - rho);
-    }
-  }
-
-  if (fCreateHisto) {
-    const Int_t Ntracks = fTracks->GetEntriesFast();
-    fHistCentrality->Fill(fCent);
-    fHistNjetvsCent->Fill(fCent, NjetAcc);
-    fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
-    PostData(1, fOutputList);
   }
-}      
 
-//________________________________________________________________________
-void AliAnalysisTaskRho::ExecOnce() 
-{
-  // Initialize some settings that need to be determined in UserExec.
-
-  AliAnalysisTaskRhoBase::ExecOnce();
-
-  if (fRhoScaled) {
-    // add rho to event if not yet there
-    if (!(InputEvent()->FindListObject(fRhoScaledName))) {
-      InputEvent()->AddObject(fRhoScaled);
-    } else {
-      AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoScaledName.Data()));
-      return;
-    }
-  }
-
-  if (fCreateHisto) {
-    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
-    if (!fTracks) {
-      AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data() ));
-      return;
-    }
-  }
-
-  fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
-  if (!fJets) {
-    AliError(Form("%s: Pointer to jets %s == 0", GetName(), fJetsName.Data() ));
-    return;
-  }
-}
-
-//________________________________________________________________________
-Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
-{
-  // Get scale factor.
-
-  Double_t scale = 1;
-  if (fScaleFunction)
-    scale = fScaleFunction->Eval(cent);
-  return scale;
-}
+  return kTRUE;
+} 
index 81182b725ecf432c254064c450986c0855408f14..b54a79786fb6f345c2dc235340590b3017bd8e1a 100644 (file)
@@ -3,14 +3,6 @@
 
 // $Id$
 
-class TClonesArray;
-class TF1;
-class TH1F;
-class TH2F;
-class TList;
-class TString;
-class AliRhoParameter;
-
 #include "AliAnalysisTaskRhoBase.h"
 
 class AliAnalysisTaskRho : public AliAnalysisTaskRhoBase {
@@ -19,62 +11,17 @@ class AliAnalysisTaskRho : public AliAnalysisTaskRhoBase {
   AliAnalysisTaskRho();
   AliAnalysisTaskRho(const char *name, Bool_t histo=kFALSE);
   virtual ~AliAnalysisTaskRho() {}
-  
-  void                   UserCreateOutputObjects();
-  void                   UserExec(Option_t*);
 
-  const char            *GetRhoScaled() const                                  { return fRhoScaledName ;                   }
-  void                   SetAreaCut(Double_t a = 0.05)                         { fAreaCut       = a    ;                   }
-  void                   SetAreaEmcCut(Double_t a = 0.99)                      { fAreaEmcCut    = a    ;                   }
-  void                   SetExcludeLeadJets(UInt_t n)                          { fNExclLeadJets = n    ;                   }
-  void                   SetJetEta(Double_t emin, Double_t emax)               { fEtaMin        = emin ; fEtaMax = emax  ; }
-  void                   SetJetPhi(Double_t pmin, Double_t pmax)               { fPhiMin        = pmin ; fPhiMax = pmax  ; }
-  void                   SetJetsName(const char *n)                            { fJetsName      = n    ;                   }
-  void                   SetRhoName(const char *name)                          { fRhoName       = name ; 
-                                                                                 fRhoScaledName = name ;
-                                                                                 fRhoScaledName += "_Scaled";              }
-  void                   SetScaleFunction(TF1* sf)                             { fScaleFunction = sf   ;                   }
-  void                   SetTracksName(const char *n)                          { fTracksName    = n    ;                   }
+  void             SetExcludeLeadJets(UInt_t n)    { fNExclLeadJets = n    ; }
 
  protected:
-  virtual void           ExecOnce();
-  virtual Double_t       GetScaleFactor(Double_t cent);
+  Bool_t           Run();
 
-  TString                fTracksName;                    // name of track collection
-  TString                fJetsName;                      // name of jet collection
-  TString                fRhoScaledName;                 // name of scaled rho object
-  Double_t               fPhiMin;                        // minimum phi
-  Double_t               fPhiMax;                        // maximum phi
-  Double_t               fEtaMin;                        // minimum eta
-  Double_t               fEtaMax;                        // maximum eta
-  Double_t               fAreaCut;                       // minimum cut on jet area
-  Double_t               fAreaEmcCut;                    // minimum cut on jet emcal area
-  UInt_t                 fNExclLeadJets;                 // number of leading jets to be excluded from the median calculation
-  TF1                   *fScaleFunction;                 // pre-computed scale factor as a function of centrality
-  Bool_t                 fCreateHisto;                   // whether or not create histograms
-  TClonesArray          *fTracks;                        //!ptr to input tracks
-  TClonesArray          *fJets;                          //!ptr to input jets
-  TList                 *fOutputList;                    //!output list
-  TH1F                  *fHistCentrality;                //!centrality distribution
-  TH1F                  *fHistJetPt;                     //!jet pt distribution
-  TH1F                  *fHistJetArea;                   //!jet area
-  TH2F                  *fHistRhovsCent;                 //!rho vs. centrality
-  TH2F                  *fHistDeltaRhovsCent;            //!delta rho vs. centrality
-  TH2F                  *fHistDeltaRhoScalevsCent;       //!delta rhoscaled vs. centrality
-  TH2F                  *fHistJetPtvsCent;               //!jet pt vs. centrality
-  TH2F                  *fHistJetAreavsCent;             //!jet area vs. centrality
-  TH2F                  *fHistNjetvsCent;                //!no. of jets vs. centrality
-  TH2F                  *fHistRhovsNtrack;               //!rho vs. no. of tracks
-  TH2F                  *fHistDeltaRhovsNtrack;          //!delta rho vs. no. of tracks
-  TH2F                  *fHistDeltaRhoScalevsNtrack;     //!delta rho scaled vs. no. of tracks
-  TH2F                  *fHistJetPtvsNtrack;             //!jet pt vs. no. of tracks
-  TH2F                  *fHistJetAreavsNtrack;           //!jet area vs. no. of tracks
-  TH2F                  *fHistNjetvsNtrack;              //!no. of jets vs. no. of tracks
-  AliRhoParameter       *fRhoScaled;                     //!per event scaled rho
+  UInt_t           fNExclLeadJets;                 // number of leading jets to be excluded from the median calculation
 
   AliAnalysisTaskRho(const AliAnalysisTaskRho&);             // not implemented
   AliAnalysisTaskRho& operator=(const AliAnalysisTaskRho&);  // not implemented
   
-  ClassDef(AliAnalysisTaskRho, 7); // Rho task
+  ClassDef(AliAnalysisTaskRho, 8); // Rho task
 };
 #endif
index 4b1837c75ec48bf85674056a25d7acc35161e60d..7ef498733e8e888290b2bb208dfbe6298b861a87 100644 (file)
 // $Id$
 //
-// Calculation of rho, method: sum of all particle pt / full acceptance area.
+// Calculation of rho, method: median all particle pt / multiplicity density.
 //
 // Authors: S. Aiola
 
 #include "AliAnalysisTaskRhoAverage.h"
 
 #include <TClonesArray.h>
-#include <TList.h>
 #include <TMath.h>
 
-#include "AliAnalysisManager.h"
-#include "AliCentrality.h"
-#include "AliEmcalJet.h"
 #include "AliLog.h"
 #include "AliRhoParameter.h"
 #include "AliVCluster.h"
-#include "AliVEventHandler.h"
 #include "AliVTrack.h"
 
 ClassImp(AliAnalysisTaskRhoAverage)
 
 //________________________________________________________________________
 AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage() : 
-  AliAnalysisTaskRhoBase(),
-  fTracksName(),
-  fClustersName(),
-  fJetsName(),
-  fEtaMin(0),
-  fEtaMax(0),
-  fPhiMin(0),
-  fPhiMax(0),
-  fPtMin(0),
-  fClusters(0),
-  fJets(0),
-  fTracks(0)
+  AliAnalysisTaskRhoBase("AliAnalysisTaskRhoAverage"),
+  fRhoType(0),
+  fNExclLeadPart(0)
 {
   // Default constructor.
 }
 
 //________________________________________________________________________
-AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage(const char *name) :
-  AliAnalysisTaskRhoBase(name),
-  fTracksName("tracks"),
-  fClustersName("caloClusters"),
-  fJetsName("KtJets"),
-  fEtaMin(-0.9),
-  fEtaMax(0.9),
-  fPhiMin(0),
-  fPhiMax(2 * TMath::Pi()),
-  fPtMin(0.15),
-  fClusters(0),
-  fJets(0),
-  fTracks(0)
+AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage(const char *name, Bool_t histo) :
+  AliAnalysisTaskRhoBase(name, histo),
+  fRhoType(0),
+  fNExclLeadPart(0)
 {
   // Constructor.
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskRhoAverage::UserExec(Option_t *
+Bool_t AliAnalysisTaskRhoAverage::Run(
 {
-  // Main loop, called for each event.
-  
-  if (!fIsInit) {
-    ExecOnce();
-    fIsInit = 1;
-  }
+  // Run the analysis.
 
-  fRho->SetVal(-1);
+  static Double_t rhovec[9999];
+  Int_t NpartAcc = 0;
 
-  Double_t rho = 0;
-  
-  Int_t Ntracks = 0;
-  if (fTracks) 
-    Ntracks = fTracks->GetEntriesFast();
-
-  Int_t Nclusters = 0;
-  if (fClusters)
-    Nclusters = fClusters->GetEntriesFast();
-
-  Int_t Njets = 0;
-  if (fJets)
-    Njets = fJets->GetEntriesFast();
-
-  Double_t maxJetPt = 0;
-  Int_t maxJetId = -1;
-  AliEmcalJet *maxJet = 0;
-  for (Int_t ij = 0; ij < Njets; ij++) {
+  Int_t   maxPartIds[] = {-1, -1};
+  Float_t maxPartPts[] = { 0,  0};
+
+  // push all jets within selected acceptance into stack
+
+  if (fNExclLeadPart > 0) {
+
+    if (fTracks && (fRhoType == 0 || fRhoType == 1)) {
       
-    AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
-  
-    if (!jet) {
-      AliError(Form("%s: Could not receive jet %d", GetName(), ij));
-      continue;
-    } 
-  
-    if (jet->Pt() > maxJetPt) {
-      maxJetPt = jet->Pt();
-      maxJetId = ij;
+      const Int_t Ntracks = fTracks->GetEntriesFast();
+      
+      for (Int_t it = 0; it < Ntracks; ++it) {
+       
+       AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
+       
+       if (!track) {
+         AliError(Form("%s: Could not receive track %d", GetName(), it));
+         continue;
+       } 
+       
+       if (!AcceptTrack(track))
+         continue;
+       
+       if (track->Pt() > maxPartPts[0]) {
+         maxPartPts[1] = maxPartPts[0];
+         maxPartIds[1] = maxPartIds[0];
+         maxPartPts[0] = track->Pt();
+         maxPartIds[0] = it+1;
+       } 
+       else if (track->Pt() > maxPartPts[1]) {
+         maxPartPts[1] = track->Pt();
+         maxPartIds[1] = it+1;
+       }
+      }
     }
-  }
 
-  if (maxJetId >= 0)
-    maxJet = static_cast<AliEmcalJet*>(fJets->At(maxJetId));
+    if (fCaloClusters && (fRhoType == 0 || fRhoType == 2)) {
 
-  for (Int_t it = 0; it < Ntracks; ++it) {
+      const Int_t Nclusters = fCaloClusters->GetEntriesFast();
       
-    AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
+      for (Int_t ic = 0; ic < Nclusters; ++ic) {
+       
+       AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
+       
+       if (!cluster) {
+         AliError(Form("%s: Could not receive cluster %d", GetName(), ic));
+         continue;
+       } 
+       
+       if (!AcceptCluster(cluster))
+         continue;
+       
+       TLorentzVector nPart;
+       cluster->GetMomentum(nPart, fVertex);
+       
+       if (nPart.Pt() > maxPartPts[0]) {
+         maxPartPts[1] = maxPartPts[0];
+         maxPartIds[1] = maxPartIds[0];
+         maxPartPts[0] = nPart.Pt();
+         maxPartIds[0] = -ic-1;
+       } 
+       else if (nPart.Pt() > maxPartPts[1]) {
+         maxPartPts[1] = nPart.Pt();
+         maxPartIds[1] = -ic-1;
+       }
+      }
+    }
+    if (fNExclLeadPart < 2) {
+      maxPartIds[1] = -1;
+      maxPartPts[1] = 0;
+    }
+  }
   
-    if (!track) {
-      AliError(Form("%s: Could not receive track %d", GetName(), it));
-      continue;
-    } 
-
-    if (track->Eta() < fEtaMin || track->Eta() > fEtaMax || track->Phi() < fPhiMin || track->Phi() > fPhiMax)
-      continue;
+  if (fTracks && (fRhoType == 0 || fRhoType == 1)) {
 
-    if (track->Pt() < fPtMin)
-      continue;
+    const Int_t Ntracks = fTracks->GetEntriesFast();
 
-    if (maxJet && IsJetTrack(maxJet, it))
-      continue;
-
-    rho += track->Pt();
-  }
-  
-  Double_t vertex[] = {0, 0, 0};
-  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+    for (Int_t it = 0; it < Ntracks; ++it) {
 
-  for (Int_t ic = 0; ic < Nclusters; ++ic) {
+      // exlcuding lead particles
+      if (it == maxPartIds[0]-1 || it == maxPartIds[1]-1)
+       continue;
       
-    AliVCluster *cluster = static_cast<AliVCluster*>(fClusters->At(ic));
+      AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
   
-    if (!cluster) {
-      AliError(Form("%s: Could not receive cluster %d", GetName(), ic));
-      continue;
-    } 
+      if (!track) {
+       AliError(Form("%s: Could not receive track %d", GetName(), it));
+       continue;
+      
 
-    Float_t pos[3];
-    cluster->GetPosition(pos);
-    TVector3 clusVec(pos);
+      if (!AcceptTrack(track))
+       continue;
 
-    if (clusVec.Eta() < fEtaMin || clusVec.Eta() > fEtaMax || 
-        clusVec.Phi() < fPhiMin || clusVec.Phi() > fPhiMax)
-      continue;
-
-    TLorentzVector nPart;
-    cluster->GetMomentum(nPart, const_cast<Double_t*>(vertex));
-
-    if (nPart.Et() < fPtMin)
-      continue;
-
-    if (maxJet && IsJetCluster(maxJet, ic))
-      continue;
-
-    rho += nPart.Et();
+      rhovec[NpartAcc] = track->Pt();
+      ++NpartAcc;
+    }
   }
-  Double_t area = (fEtaMax - fEtaMin) * (fPhiMax - fPhiMin);
 
-  if (maxJet)
-    area -= maxJet->Area();
+  if (fCaloClusters && (fRhoType == 0 || fRhoType == 2)) {
 
-  if (area>0) {
-    rho /= area;
-    fRho->SetVal(rho);
-  } else {
-    AliError(Form("%s: Area negative %f", GetName(), area));
-  }
-}      
+    const Int_t Nclusters = fCaloClusters->GetEntriesFast();
 
-//________________________________________________________________________
-void AliAnalysisTaskRhoAverage::ExecOnce() 
-{
-  // Initialize some settings that need to be determined in UserExec.
+    for (Int_t ic = 0; ic < Nclusters; ++ic) {
 
-  AliAnalysisTaskRhoBase::ExecOnce();
+      // exlcuding lead particles
+      if (ic == -maxPartIds[0]-1 || ic == -maxPartIds[1]-1)
+       continue;
+      
+      AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
+      
+      if (!cluster) {
+       AliError(Form("%s: Could not receive cluster %d", GetName(), ic));
+       continue;
+      } 
+      
+      if (!AcceptCluster(cluster))
+       continue;
+      
+      TLorentzVector nPart;
+      cluster->GetMomentum(nPart, fVertex);
 
-  if (!fClustersName.IsNull()) {
-    fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fClustersName));
-    if (!fClusters) {
-      AliError(Form("%s: Pointer to jets %s == 0", GetName(), fClustersName.Data() ));
-      return;
+      rhovec[NpartAcc] = nPart.Pt();
+      ++NpartAcc;
     }
   }
 
-  if (!fTracksName.IsNull()) {
-    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
-    if (!fTracks) {
-      AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data() ));
-      return;
-    }
-  }
+  Double_t rho = TMath::Median(NpartAcc, rhovec);
+  Double_t area = (fTrackMaxEta - fTrackMinEta) * (fTrackMaxPhi - fTrackMinPhi);
 
-  if (!fJetsName.IsNull()) {
-    fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
-    if (!fJets) {
-      AliError(Form("%s: Pointer to jets %s == 0", GetName(), fJetsName.Data() ));
-      return;
-    }
+  if (area > 0) {
+    rho *= NpartAcc / area;
+    fRho->SetVal(rho);
+  } 
+  else {
+    AliError(Form("%s: Area <= 0 %f", GetName(), area));
+    return kFALSE;
   }
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskRhoAverage::IsJetTrack(AliEmcalJet* jet, Int_t itrack) const
-{
-  // Return true if track is in jet.
 
-  for (Int_t i = 0; i < jet->GetNumberOfTracks(); i++) {
-    Int_t ijettrack = jet->TrackAt(i);
-    if (ijettrack == itrack)
-      return kTRUE;
+  if (fScaleFunction) {
+    Double_t rhoScaled = rho * GetScaleFactor(fCent);
+    fRhoScaled->SetVal(rhoScaled);
   }
-  return kFALSE;
-}
 
-//________________________________________________________________________
-Bool_t AliAnalysisTaskRhoAverage::IsJetCluster(AliEmcalJet* jet, Int_t iclus) const
-{
-  // Return true if cluster is in jet.
-
-  for (Int_t i = 0; i < jet->GetNumberOfClusters(); i++) {
-    Int_t ijetclus = jet->ClusterAt(i);
-    if (ijetclus == iclus)
-      return kTRUE;
-  }
-  return kFALSE;
+  return kTRUE;
 }
index bbf91560f0af1fa559873f2dad004de2e304336a..e85682ce413da845602e21a988af56667ce228b5 100644 (file)
@@ -3,48 +3,27 @@
 
 // $Id$
 
-class TClonesArray;
-class TList;
-class AliEmcalJet;
-
 #include "AliAnalysisTaskRhoBase.h"
 
 class AliAnalysisTaskRhoAverage : public AliAnalysisTaskRhoBase {
 
  public:
   AliAnalysisTaskRhoAverage();
-  AliAnalysisTaskRhoAverage(const char *name);
+  AliAnalysisTaskRhoAverage(const char *name, Bool_t histo=kFALSE);
   virtual ~AliAnalysisTaskRhoAverage() {}
-  
-  void                   UserExec(Option_t*);
 
-  void                   SetClustersName(const char *n)                        { fClustersName  = n    ; }
-  void                   SetEtaLimits(Double_t emin, Double_t emax)            { fEtaMin        = emin ; fEtaMax = emax  ; }
-  void                   SetJetsName(const char *n)                            { fJetsName      = n    ; }  
-  void                   SetPhiLimits(Double_t pmin, Double_t pmax)            { fPhiMin        = pmin ; fPhiMax = pmax  ; }
-  void                   SetPtMin(Double_t pt)                                 { fPtMin         = pt   ; }
-  void                   SetTracksName(const char *n)                          { fTracksName    = n    ; }
+  void             SetRhoType(Int_t t)             { fRhoType       = t    ; }
+  void             SetExcludeLeadPart(UInt_t n)    { fNExclLeadPart = n    ; }
   
  protected:
-  void                   ExecOnce();
-  Bool_t                 IsJetCluster(AliEmcalJet* jet, Int_t iclus) const;
-  Bool_t                 IsJetTrack(AliEmcalJet* jet, Int_t itrack)  const;
+  Bool_t           Run();
 
-  TString                fTracksName;                    // name of track collection
-  TString                fClustersName;                  // name of clusters collection
-  TString                fJetsName;                      // name of jet collection
-  Double_t               fEtaMin;                        // minimum eta
-  Double_t               fEtaMax;                        // maximum eta
-  Double_t               fPhiMin;                        // minimum phi
-  Double_t               fPhiMax;                        // maximum phi
-  Double_t               fPtMin;                         // minimum pt
-  TClonesArray          *fClusters;                      //!input clusters
-  TClonesArray          *fJets;                          //!input jets
-  TClonesArray          *fTracks;                        //!input tracks
+  Int_t            fRhoType       ;// rho type: 0 = charged+neutral, 1 = charged, 2 = neutral
+  UInt_t           fNExclLeadPart ;// number of leading particles to be excluded from the median calculation
 
   AliAnalysisTaskRhoAverage(const AliAnalysisTaskRhoAverage&);             // not implemented
   AliAnalysisTaskRhoAverage& operator=(const AliAnalysisTaskRhoAverage&);  // not implemented
   
-  ClassDef(AliAnalysisTaskRhoAverage, 2); // Rho task, method: sum of all particle pt / full acceptance area
+  ClassDef(AliAnalysisTaskRhoAverage, 3); // Rho task
 };
 #endif
index 70be2f7f979601a46c7ef7b48b86d9b201adc91e..a3afeb6d735d311f8ee43e5f69e845f4e00a710e 100644 (file)
@@ -6,15 +6,14 @@
 // Author: S.Aiola
 
 #include <TF1.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TClonesArray.h>
 
 #include "AliAnalysisManager.h"
-#include "AliCentrality.h"
-#include "AliESDEvent.h"
-#include "AliEmcalJet.h"
 #include "AliLog.h"
 #include "AliRhoParameter.h"
-#include "AliVCluster.h"
-#include "AliVEventHandler.h"
+#include "AliEmcalJet.h"
 
 #include "AliAnalysisTaskRhoBase.h"
 
@@ -22,88 +21,247 @@ ClassImp(AliAnalysisTaskRhoBase)
 
 //________________________________________________________________________
 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase() : 
-  AliAnalysisTaskSE(),
-  fRhoName("Rho"),
-  fRhoFunction(0x0),
-  fCent(-1),
-  fRho(0),
-  fDoCent(0),
-  fIsInit(0)
+  AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoBase", kFALSE),
+  fRhoScaledName(),
+  fCompareRhoName(),
+  fCompareRhoScaledName(),
+  fRhoFunction(0),
+  fScaleFunction(0),
+  fRhoScaled(0),
+  fCompareRho(0),
+  fCompareRhoScaled(0),
+  fHistJetPtvsCent(0),
+  fHistJetAreavsCent(0),
+  fHistNjetvsCent(0),
+  fHistJetPtvsNtrack(0),
+  fHistJetAreavsNtrack(0),
+  fHistNjetvsNtrack(0),
+  fHistRhovsCent(0),
+  fHistRhoScaledvsCent(0),
+  fHistDeltaRhovsCent(0),
+  fHistDeltaRhoScalevsCent(0),
+  fHistRhovsNtrack(0),
+  fHistRhoScaledvsNtrack(0),
+  fHistDeltaRhovsNtrack(0),
+  fHistDeltaRhoScalevsNtrack(0),
+  fHistRhovsNcluster(0),
+  fHistRhoScaledvsNcluster(0)
 {
   // Constructor.
 }
 
 //________________________________________________________________________
-AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase(const char *name) :
-  AliAnalysisTaskSE(name),
-  fRhoName("Rho"),
-  fRhoFunction(0x0),
-  fCent(-1),
-  fRho(0),
-  fDoCent(0),
-  fIsInit(0)
+AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase(const char *name, Bool_t histo) :
+  AliAnalysisTaskEmcalJet(name, histo),
+  fRhoScaledName(),
+  fCompareRhoName(),
+  fCompareRhoScaledName(),
+  fRhoFunction(0),
+  fScaleFunction(0),
+  fRhoScaled(0),
+  fCompareRho(0),
+  fCompareRhoScaled(0),
+  fHistJetPtvsCent(0),
+  fHistJetAreavsCent(0),
+  fHistNjetvsCent(0),
+  fHistJetPtvsNtrack(0),
+  fHistJetAreavsNtrack(0),
+  fHistNjetvsNtrack(0),
+  fHistRhovsCent(0),
+  fHistRhoScaledvsCent(0),
+  fHistDeltaRhovsCent(0),
+  fHistDeltaRhoScalevsCent(0),
+  fHistRhovsNtrack(0),
+  fHistRhoScaledvsNtrack(0),
+  fHistDeltaRhovsNtrack(0),
+  fHistDeltaRhoScalevsNtrack(0),
+  fHistRhovsNcluster(0),
+  fHistRhoScaledvsNcluster(0)
 {
   // Constructor.
+
+  SetMakeGeneralHistograms(histo);
 }
 
 //________________________________________________________________________
 void AliAnalysisTaskRhoBase::UserCreateOutputObjects()
 {
-  // Run at beginning of task.
+  // User create output objects, called at the beginning of the analysis.
 
-  AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
-  if (!handler) {
-    AliError("Input handler not available!");
+  if (!fCreateHisto)
     return;
+
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+  fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt*2);
+  fOutput->Add(fHistRhovsCent);
+
+  if (!fTracksName.IsNull()) {
+    fHistRhovsNtrack = new TH2F("RhovsNtrack", "RhovsNtrack", 125, 0, 4000, fNbins, fMinBinPt, fMaxBinPt*2);
+    fOutput->Add(fHistRhovsNtrack);
   }
 
-  fRho = new AliRhoParameter(fRhoName, 0);
-}
+  if (!fCaloName.IsNull()) {
+    fHistRhovsNcluster = new TH2F("RhovsNcluster", "RhovsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
+    fOutput->Add(fHistRhovsNcluster);
+  }
 
-//________________________________________________________________________
-void AliAnalysisTaskRhoBase::UserExec(Option_t *) 
-{
-  // Main loop, called for each event.
+  if (!fJetsName.IsNull()) {
+    fHistJetPtvsCent            = new TH2F("JetPtvsCent",           "JetPtvsCent",           101, -1,  100,   fNbins, fMinBinPt, fMaxBinPt);
+    fHistJetAreavsCent          = new TH2F("JetAreavsCent",         "JetAreavsCent",         101, -1,  100,   30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
+    fHistNjetvsCent             = new TH2F("NjetvsCent",            "NjetvsCent",            101, -1,  100,   150, -0.5, 149.5);
+
+    fOutput->Add(fHistJetPtvsCent);
+    fOutput->Add(fHistJetAreavsCent);
+    fOutput->Add(fHistNjetvsCent);
+
+    if (!fTracksName.IsNull()) {
+      fHistJetPtvsNtrack        = new TH2F("JetPtvsNtrack",         "JetPtvsNtrack",         125,  0,  4000,  fNbins, fMinBinPt, fMaxBinPt);
+      fHistJetAreavsNtrack      = new TH2F("JetAreavsNtrack",       "JetAreavsNtrack",       125,  0,  4000,  30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
+      fHistNjetvsNtrack         = new TH2F("NjetvsNtrack",          "rNjetvsNtrack",         125,  0,  4000,  150, -0.5, 149.5);
 
-  if (!fIsInit) {
-    ExecOnce();
-    fIsInit = 1;
+      fOutput->Add(fHistJetPtvsNtrack);
+      fOutput->Add(fHistJetAreavsNtrack);
+      fOutput->Add(fHistNjetvsNtrack);
+    }
+  }
+  
+  if (!fCompareRhoName.IsNull()) {
+    fHistDeltaRhovsCent = new TH2F("DeltaRhovsCent", "DetlaRhovsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
+    fOutput->Add(fHistDeltaRhovsCent);
+    if (!fTracksName.IsNull()) {
+      fHistDeltaRhovsNtrack = new TH2F("DeltaRhovsNtrack", "DeltaRhovsNtrack", 125, 0, 4000, fNbins, -fMaxBinPt, fMaxBinPt);
+      fOutput->Add(fHistDeltaRhovsNtrack);
+    }
   }
 
-  DetermineCent();
+  if (fScaleFunction) {
+    fHistRhoScaledvsCent = new TH2F("RhoScaledvsCent", "RhoScalevsCent", 101, -1, 100, fNbins, fMinBinPt , fMaxBinPt*2);
+    fOutput->Add(fHistRhoScaledvsCent);
+
+    if (!fTracksName.IsNull()) {
+      fHistRhoScaledvsNtrack = new TH2F("RhoScaledvsNtrack", "RhoScaledvsNtrack", 125, 0, 4000, fNbins, fMinBinPt, fMaxBinPt*2);
+      fOutput->Add(fHistRhoScaledvsNtrack);
+    }
+
+    if (!fCaloName.IsNull()) {
+      fHistRhoScaledvsNcluster = new TH2F("RhoScaledvsNcluster", "RhoScaledvsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
+      fOutput->Add(fHistRhoScaledvsNcluster);
+    }
+
+    if (!fCompareRhoScaledName.IsNull()) {
+      fHistDeltaRhoScalevsCent = new TH2F("DeltaRhoScalevsCent", "DeltaRhoScalevsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
+      fOutput->Add(fHistDeltaRhoScalevsCent);
+      
+      if (!fTracksName.IsNull()) {
+       fHistDeltaRhoScalevsNtrack = new TH2F("DeltaRhoScalevsNtrack", "DeltaRhoScalevsNtrack", 125, 0, 4000, fNbins, -fMaxBinPt, fMaxBinPt);
+       fOutput->Add(fHistDeltaRhoScalevsNtrack);
+      }
+    }
+  }
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskRhoBase::Run() 
+{
+  // Run the analysis.
 
   Double_t rho = GetRhoFactor(fCent);
   fRho->SetVal(rho);
-}      
+
+  if (fScaleFunction) {
+    Double_t rhoScaled = rho * GetScaleFactor(fCent);
+    fRhoScaled->SetVal(rhoScaled);
+  }
+
+  return kTRUE;
+}
 
 //________________________________________________________________________
-void AliAnalysisTaskRhoBase::DetermineCent() 
+Bool_t AliAnalysisTaskRhoBase::FillHistograms() 
 {
-  // Determine centrality.
+  // Fill histograms.
 
-  fCent = 99; 
-  
-  if (fDoCent) {
-    AliCentrality *centrality = InputEvent()->GetCentrality();
+  Int_t Njets     = 0;
+  Int_t Ntracks   = 0;
+  Int_t Nclusters = 0;
+
+  if (fJets)
+    Njets     = fJets->GetEntries();
+  if (fTracks)
+    Ntracks   = fTracks->GetEntriesFast();
+  if (fCaloClusters)
+    Nclusters = fCaloClusters->GetEntriesFast();
+
+  Int_t NjetAcc = 0;
+
+  if (fJets) {
+
+    for (Int_t i = 0; i < Njets; ++i) {
+      
+      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(i));
+      if (!jet) {
+       AliError(Form("%s: Could not receive jet %d", GetName(), i));
+       continue;
+      } 
+      
+      if (!AcceptJet(jet))
+       continue;
+      
+      fHistJetPtvsCent->Fill(fCent, jet->Pt());
+      fHistJetAreavsCent->Fill(fCent, jet->Area());
+      
+      if (fTracks) {
+       fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
+       fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
+      }
+      
+      NjetAcc++;
+    }
     
-    if (centrality)
-      fCent = centrality->GetCentralityPercentile("V0M");
-    else
-      fCent = 99; // probably pp data
     
-    if (fCent < 0) {
-      AliWarning(Form("%s: Centrality negative: %f, assuming 99", GetName(), fCent));
-      fCent = 99;
+    fHistNjetvsCent->Fill(fCent, NjetAcc);
+    if (fTracks)
+      fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
+  }
+  
+  fHistRhovsCent->Fill(fCent, fRho->GetVal());
+
+  if (fTracks)
+    fHistRhovsNtrack->Fill(Ntracks, fRho->GetVal());
+  if (fCaloClusters)
+    fHistRhovsNcluster->Fill(Nclusters, fRho->GetVal());
+  if (fCompareRho) {
+    fHistDeltaRhovsCent->Fill(fCent, fRho->GetVal() - fCompareRho->GetVal());
+    if (fTracks)
+      fHistDeltaRhovsNtrack->Fill(Ntracks, fRho->GetVal() - fCompareRho->GetVal());
+  }
+
+  if (fScaleFunction) {
+    fHistRhoScaledvsCent->Fill(fCent, fRhoScaled->GetVal());
+    if (fTracks)
+      fHistRhoScaledvsNtrack->Fill(Ntracks, fRhoScaled->GetVal());
+    if (fCaloClusters)
+      fHistRhoScaledvsNcluster->Fill(Nclusters,  fRhoScaled->GetVal());
+    if (fCompareRhoScaled) {
+      fHistDeltaRhoScalevsCent->Fill(fCent, fRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
+      if (fTracks)
+       fHistDeltaRhoScalevsNtrack->Fill(Ntracks, fRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
     }
   }
-}
+
+  return kTRUE;
+}      
+
 
 //________________________________________________________________________
 void AliAnalysisTaskRhoBase::ExecOnce() 
 {
-  // Initialize some settings that need to be determined in UserExec.
+  // Init the analysis.
 
   // add rho to event if not yet there
+  fRho = new AliRhoParameter(fRhoName, 0);
+
   if (!(InputEvent()->FindListObject(fRhoName))) {
     InputEvent()->AddObject(fRho);
   } else {
@@ -111,42 +269,31 @@ void AliAnalysisTaskRhoBase::ExecOnce()
     return;
   }
 
-  // determine if centrality should be used
-  TString bt(GetBeamType());
-  if (bt == "A-A")
-    fDoCent = 1;
-  else fDoCent = 0;
-}
-
-//_____________________________________________________
-TString AliAnalysisTaskRhoBase::GetBeamType()
-{
-  // Get beam type : pp-AA-pA
-  // ESDs have it directly, AODs get it from hardcoded run number ranges
-  
-  AliVEvent *event = InputEvent();
-  if (!event) { 
-    AliError(Form("%s: Couldn't retrieve event!", GetName()));
-    return "";
+  if (fScaleFunction) {
+    fRhoScaled = new AliRhoParameter(fRhoScaledName, 0);
+    if (!(InputEvent()->FindListObject(fRhoScaledName))) {
+      InputEvent()->AddObject(fRhoScaled);
+    } else {
+      AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoScaledName.Data()));
+      return;
+    }
   }
 
-  TString beamType;
+  if (!fCompareRhoName.IsNull() && !fCompareRho) {
+    fCompareRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoName));
+    if (!fCompareRho) {
+      AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoName.Data()));
+    }
+  }
 
-  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
-  if (esd) {
-    const AliESDRun *run = esd->GetESDRun();
-    beamType = run->GetBeamType();
-  } else {
-    Int_t runNumber = event->GetRunNumber();
-    if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
-       (runNumber >= 166529 && runNumber <= 170593)) {  // LHC11h
-      beamType = "A-A";
-    } else {
-      beamType = "p-p";
+  if (!fCompareRhoScaledName.IsNull() && !fCompareRhoScaled) {
+    fCompareRhoScaled = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoScaledName));
+    if (!fCompareRhoScaled) {
+      AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoScaledName.Data()));
     }
   }
 
-  return beamType;    
+  AliAnalysisTaskEmcalJet::ExecOnce();
 }
 
 //________________________________________________________________________
@@ -154,8 +301,19 @@ Double_t AliAnalysisTaskRhoBase::GetRhoFactor(Double_t cent)
 {
   // Return rho per centrality.
 
-  Double_t rho = -1;
+  Double_t rho = 0;
   if (fRhoFunction)
     rho = fRhoFunction->Eval(cent);
   return rho;
 }
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskRhoBase::GetScaleFactor(Double_t cent)
+{
+  // Get scale factor.
+
+  Double_t scale = 1;
+  if (fScaleFunction)
+    scale = fScaleFunction->Eval(cent);
+  return scale;
+}
index 71bf4abd057878b968016e025a32f9681f39e0a5..9b20a0b08c269020063c8bf93d24b3dae0205719 100644 (file)
@@ -5,39 +5,68 @@
 
 class TString;
 class TF1;
+class TH1F;
+class TH2F;
 class AliRhoParameter;
 
-#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskEmcalJet.h"
 
-class AliAnalysisTaskRhoBase : public AliAnalysisTaskSE {
+class AliAnalysisTaskRhoBase : public AliAnalysisTaskEmcalJet {
  public:
   AliAnalysisTaskRhoBase();
-  AliAnalysisTaskRhoBase(const char *name);
+  AliAnalysisTaskRhoBase(const char *name, Bool_t histo=kFALSE);
   virtual ~AliAnalysisTaskRhoBase() {}
-  
+
   void                   UserCreateOutputObjects();
-  void                   UserExec(Option_t*);
 
-  const char            *GetRhoName() const                                    { return fRhoName       ; }
-  void                   SetRhoFunction(TF1* rf)                               { fRhoFunction   = rf   ; }
-  void                   SetRhoName(const char *name)                          { fRhoName       = name ; }
+  void                   SetRhoName(const char *name)                          { fRhoName              = name ; 
+                                                                                 fRhoScaledName        = Form("%s_Scaled",name) ; }
+  void                   SetCompareRhoName(const char *name)                   { fCompareRhoName       = name ;                   }
+  void                   SetCompareRhoScaledName(const char *name)             { fCompareRhoScaledName = name ;                   }
+  void                   SetScaleFunction(TF1* sf)                             { fScaleFunction        = sf   ;                   }
+  void                   SetRhoFunction(TF1* rf)                               { fRhoFunction          = rf   ;                   }
 
  protected:
-  virtual void           DetermineCent();
-  virtual void           ExecOnce();
-  TString                GetBeamType();
+  void                   ExecOnce();
+  Bool_t                 Run();
+  Bool_t                 FillHistograms();
+
   virtual Double_t       GetRhoFactor(Double_t cent);
+  virtual Double_t       GetScaleFactor(Double_t cent);
 
-  TString                fRhoName;                       // name of rho
+  TString                fRhoScaledName;                 // name of scaled rho object
+  TString                fCompareRhoName;                // name of rho object to compare
+  TString                fCompareRhoScaledName;          // name of scaled rho object to compare
   TF1                   *fRhoFunction;                   // pre-computed rho as a function of centrality
-  Double_t               fCent;                          //!event centrality
-  AliRhoParameter       *fRho;                           //!per event calculated rho
-  Bool_t                 fDoCent;                        //!==1 then do centrality
-  Bool_t                 fIsInit;                        //!==1 then do init
+  TF1                   *fScaleFunction;                 // pre-computed scale factor as a function of centrality
+
+  AliRhoParameter       *fRhoScaled;                     //!scaled rho object
+  AliRhoParameter       *fCompareRho;                    //!rho object to compare
+  AliRhoParameter       *fCompareRhoScaled;              //!scaled rho object to compare
+
+  TH2F                  *fHistJetPtvsCent;               //!jet pt vs. centrality
+  TH2F                  *fHistJetAreavsCent;             //!jet area vs. centrality
+  TH2F                  *fHistNjetvsCent;                //!no. of jets vs. centrality
+  TH2F                  *fHistJetPtvsNtrack;             //!jet pt vs. no. of tracks
+  TH2F                  *fHistJetAreavsNtrack;           //!jet area vs. no. of tracks
+  TH2F                  *fHistNjetvsNtrack;              //!no. of jets vs. no. of tracks
+
+  TH2F                  *fHistRhovsCent;                 //!rho vs. centrality
+  TH2F                  *fHistRhoScaledvsCent;           //!rhoscaled vs. centrality
+  TH2F                  *fHistDeltaRhovsCent;            //!delta rho vs. centrality
+  TH2F                  *fHistDeltaRhoScalevsCent;       //!delta rhoscaled vs. centrality
+
+  TH2F                  *fHistRhovsNtrack;               //!rho vs. no. of tracks
+  TH2F                  *fHistRhoScaledvsNtrack;         //!rhoscaled vs. no. of tracks
+  TH2F                  *fHistDeltaRhovsNtrack;          //!delta rho vs. no. of tracks
+  TH2F                  *fHistDeltaRhoScalevsNtrack;     //!delta rho scaled vs. no. of tracks
+  TH2F                  *fHistRhovsNcluster;             //!rho vs. no. of clusters
+  TH2F                  *fHistRhoScaledvsNcluster;       //!rhoscaled vs. no. of clusters
 
   AliAnalysisTaskRhoBase(const AliAnalysisTaskRhoBase&);             // not implemented
   AliAnalysisTaskRhoBase& operator=(const AliAnalysisTaskRhoBase&);  // not implemented
   
-  ClassDef(AliAnalysisTaskRhoBase, 3); // Rho base task
+  ClassDef(AliAnalysisTaskRhoBase, 4); // Rho base task
 };
 #endif
index 8f6798ee06f8f6bac22189f8bab91ebd04f3fba0..adb74c44a8e55d5771f15c6d86be1f4e92822871 100644 (file)
@@ -45,7 +45,6 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fNTrials(0),
   fMCTracks(0),
   fMCJets(0),
-  fHistZVertex(0),
   fHistNTrials(0),
   fHistEvents(0),
   fHistMCJetPhiEta(0),
@@ -72,6 +71,8 @@ AliJetResponseMaker::AliJetResponseMaker() :
   for (Int_t i = 0; i < 11; i++) {
     fHistEventWeight[i] = 0;
   }
+
+  SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
@@ -95,7 +96,6 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fNTrials(0),
   fMCTracks(0),
   fMCJets(0),
-  fHistZVertex(0),
   fHistNTrials(0),
   fHistEvents(0),
   fHistMCJetPhiEta(0),
@@ -122,6 +122,8 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   for (Int_t i = 0; i < 11; i++) {
     fHistEventWeight[i] = 0;
   }
+
+  SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
@@ -135,18 +137,11 @@ void AliJetResponseMaker::UserCreateOutputObjects()
 {
   // Create user objects.
 
-  OpenFile(1);
-  fOutput = new TList();
-  fOutput->SetOwner();
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 
   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 
-  fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
-  fHistZVertex->GetXaxis()->SetTitle("z");
-  fHistZVertex->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistZVertex);
-
   fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
   fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
   fHistNTrials->GetYaxis()->SetTitle("trials");
@@ -234,7 +229,7 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
   fOutput->Add(fHistClosestDeltaPhi);
 
-  fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fMaxEta - fMinEta) * 20, fMinEta * 2, fMaxEta * 2);
+  fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2);
   fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
   fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
   fOutput->Add(fHistClosestDeltaEta);
@@ -268,7 +263,7 @@ void AliJetResponseMaker::UserCreateOutputObjects()
 }
 
 //________________________________________________________________________
-Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet, Bool_t /*bias*/, Bool_t /*upCut*/) const
+Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
 {   
   // Return true if jet is accepted.
 
@@ -317,16 +312,16 @@ void AliJetResponseMaker::ExecOnce()
   AliAnalysisTaskEmcalJet::ExecOnce();
 
   if (fMCFiducial) {
-    fMCminEta = fMinEta;
-    fMCmaxEta = fMaxEta;
-    fMCminPhi = fMinPhi;
-    fMCmaxPhi = fMaxPhi;
+    fMCminEta = fJetMinEta;
+    fMCmaxEta = fJetMaxEta;
+    fMCminPhi = fJetMinPhi;
+    fMCmaxPhi = fJetMaxPhi;
   }
   else {
-    fMCminEta = fMinEta - fJetRadius;
-    fMCmaxEta = fMaxEta + fJetRadius;
-    fMCminPhi = fMinPhi - fJetRadius;
-    fMCmaxPhi = fMaxPhi + fJetRadius;
+    fMCminEta = fJetMinEta - fJetRadius;
+    fMCmaxEta = fJetMaxEta + fJetRadius;
+    fMCminPhi = fJetMinPhi - fJetRadius;
+    fMCmaxPhi = fJetMaxPhi + fJetRadius;
   }
 }
 
@@ -436,7 +431,7 @@ void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bo
       continue;
 
     if (!mc) {
-      if (jet1->Eta() < fMinEta || jet1->Eta() > fMaxEta || jet1->Phi() < fMinPhi || jet1->Phi() > fMaxPhi)
+      if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
        continue;
     }
     else {
@@ -457,7 +452,7 @@ void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bo
        continue;
 
       if (mc) {
-       if (jet2->Eta() < fMinEta || jet2->Eta() > fMaxEta || jet2->Phi() < fMinPhi || jet2->Phi() > fMaxPhi)
+       if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
          continue;
       }
       else {
@@ -489,7 +484,6 @@ Bool_t AliJetResponseMaker::FillHistograms()
   fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
   if (fEventWeightHist)
     fHistEventWeight[fPtHardBin]->Fill(fPythiaHeader->EventWeight());
-  fHistZVertex->Fill(fVertex[2]);
 
   const Int_t nMCJets = fMCJets->GetEntriesFast();
 
@@ -552,7 +546,7 @@ Bool_t AliJetResponseMaker::FillHistograms()
 
     if (!AcceptBiasJet(jet))
       continue;
-    if (jet->Eta() < fMinEta || jet->Eta() > fMaxEta || jet->Phi() < fMinPhi || jet->Phi() > fMaxPhi)
+    if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi)
       continue;
     
     fHistMCJetsPtFiducial->Fill(jet->Pt(), fEventWeight);
@@ -576,7 +570,7 @@ Bool_t AliJetResponseMaker::FillHistograms()
       continue;
     if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
       continue;
-    if (jet->Eta() < fMinEta || jet->Eta() > fMaxEta || jet->Phi() < fMinPhi || jet->Phi() > fMaxPhi)
+    if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi)
       continue;
 
     if (!jet->MatchedJet()) {
index d57c711368cb9d16f78b6a52141535a713c342ac..b25f0ceae5fa0f3698bc4d67a9eb197d3a81bfa0 100644 (file)
@@ -29,7 +29,7 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
 
  protected:
   Bool_t                      IsEventSelected();
-  Bool_t                      AcceptJet(AliEmcalJet* jet, Bool_t /*bias*/ = kFALSE, Bool_t /*upCut*/ = kFALSE)   const;
+  Bool_t                      AcceptJet(AliEmcalJet* jet)   const;
   void                        ExecOnce();
   void                        DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t mc);
   Bool_t                      FillHistograms();
@@ -56,7 +56,6 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   TClonesArray               *fMCTracks;                  //!MC particles
   TClonesArray               *fMCJets;                    //!MC jets
   // General histograms
-  TH1F                       *fHistZVertex;               //!Z vertex position
   TH1F                       *fHistNTrials;               //!total number of trials per pt hard bin
   TH1F                       *fHistEvents;                //!total number of events per pt hard bin
   TH1F                       *fHistEventWeight[11];       //!event weight
@@ -86,6 +85,6 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   AliJetResponseMaker(const AliJetResponseMaker&);            // not implemented
   AliJetResponseMaker &operator=(const AliJetResponseMaker&); // not implemented
 
-  ClassDef(AliJetResponseMaker, 7) // Jet response matrix producing task
+  ClassDef(AliJetResponseMaker, 8) // Jet response matrix producing task
 };
 #endif
index 69309722470b51af5febbb38f160662bd6564e3e..e221bc2324d5a711274208c34046a49d7f7ca1b7 100644 (file)
@@ -9,6 +9,7 @@
 #include <TClonesArray.h>
 #include <TH1F.h>
 #include <TH2F.h>
+#include <TH3F.h>
 #include <TList.h>
 #include <TLorentzVector.h>
 #include <TRandom3.h>
@@ -31,146 +32,54 @@ ClassImp(AliAnalysisTaskSAJF)
 //________________________________________________________________________
 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() : 
   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF", kTRUE),
-  fMCAna(kFALSE),
-  fMinRC2LJ(-1),
-  fEmbJetsName("EmbJets"),
-  fEmbTracksName(""),
-  fEmbCaloName(""),
-  fRandTracksName("TracksRandomized"),
-  fRandCaloName("CaloClustersRandomized"),
-  fRhoName("Rho"),
-  fRCperEvent(-1),
-  fEmbJets(0),
-  fEmbTracks(0),
-  fEmbCaloClusters(0),
-  fRandTracks(0),
-  fRandCaloClusters(0),
-  fRho(0),
-  fRhoVal(0),
-  fEmbeddedClusterId(-1),
-  fEmbeddedTrackId(-1),
-  fHistCentrality(0),
-  fHistDeltaVectorPt(0),
-  fHistRhoVSleadJetPt(0),
-  fHistRCPhiEta(0),
-  fHistRCPtExLJVSDPhiLJ(0),
-  fHistRhoVSRCPt(0),
-  fHistEmbJetPhiEta(0),
-  fHistEmbPartPhiEta(0),
-  fHistRhoVSEmbBkg(0)
+  fLeadingHadronType(0),
+  fHistRhoVSleadJetPt(0)
 
 {
   // Default constructor.
 
   for (Int_t i = 0; i < 4; i++) {
     fHistEvents[i] = 0;
-    fHistTracksPt[i] = 0;
-    fHistClustersPt[i] = 0;
-    fHistJetPhiEta[i] = 0;
-    fHistJetsPt[i] = 0;
-    fHistJetsPtArea[i] = 0;
     fHistLeadingJetPt[i] = 0;
     fHist2LeadingJetPt[i] = 0;
+    fHistLeadingJetCorrPt[i] = 0;
+    fHistJetPhiEta[i] = 0;
+    fHistJetsPtArea[i] = 0;
+    fHistJetsCorrPtArea[i] = 0;
     fHistJetsNEFvsPt[i] = 0;
     fHistJetsZvsPt[i] = 0;
-    fHistMaxTrackPtvsJetPt[i] = 0;
-    fHistMaxClusPtvsJetPt[i] = 0;
-    fHistMaxPartPtvsJetPt[i] = 0;
-    fHistMaxTrackPtvsJetCorrPt[i] = 0;
-    fHistMaxClusPtvsJetCorrPt[i] = 0;
-    fHistMaxPartPtvsJetCorrPt[i] = 0;
     fHistConstituents[i] = 0;
-    fHistRho[i] = 0;
-    fHistJetsCorrPt[i] = 0;
-    fHistJetsCorrPtArea[i] = 0;
-    fHistLeadingJetCorrPt[i] = 0;
-    fHistRCPtRigid[i] = 0;
-    fHistRCPt[i] = 0;
-    fHistRCPtExLJ[i] = 0;
-    fHistRCPtRand[i] = 0;
-    fHistDeltaPtRCRigid[i] = 0;
-    fHistDeltaPtRC[i] = 0;
-    fHistDeltaPtRCExLJ[i] = 0;
-    fHistDeltaPtRCRand[i] = 0;
-    fHistEmbNotFoundPhiEta[i] = 0;
-    fHistEmbJetsPt[i] = 0;
-    fHistEmbJetsCorrPt[i] = 0;
-    fHistEmbPartPt[i] = 0;
-    fHistDistEmbPartJetAxis[i] = 0;
-    fHistDeltaPtEmb[i] = 0;
+    fHistTracksJetPt[i] = 0;
+    fHistClustersJetPt[i] = 0;
   }
+
+  SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) : 
   AliAnalysisTaskEmcalJet(name, kTRUE),
-  fMCAna(kFALSE),
-  fMinRC2LJ(-1),
-  fEmbJetsName("EmbJets"),
-  fEmbTracksName(""),
-  fEmbCaloName(""),
-  fRandTracksName("TracksRandomized"),
-  fRandCaloName("CaloClustersRandomized"),
-  fRhoName("Rho"), 
-  fRCperEvent(-1),
-  fEmbJets(0),
-  fEmbTracks(0),
-  fEmbCaloClusters(0),
-  fRandTracks(0),
-  fRandCaloClusters(0),
-  fRho(0),
-  fRhoVal(0),
-  fEmbeddedClusterId(-1),
-  fEmbeddedTrackId(-1),
-  fHistCentrality(0),
-  fHistDeltaVectorPt(0),
-  fHistRhoVSleadJetPt(0),
-  fHistRCPhiEta(0),
-  fHistRCPtExLJVSDPhiLJ(0),
-  fHistRhoVSRCPt(0),
-  fHistEmbJetPhiEta(0),
-  fHistEmbPartPhiEta(0),
-  fHistRhoVSEmbBkg(0)
+  fLeadingHadronType(0),
+  fHistRhoVSleadJetPt(0)
 {
   // Standard constructor.
 
   for (Int_t i = 0; i < 4; i++) {
     fHistEvents[i] = 0;
-    fHistTracksPt[i] = 0;
-    fHistClustersPt[i] = 0;
-    fHistJetPhiEta[i] = 0;
-    fHistJetsPt[i] = 0;
-    fHistJetsPtArea[i] = 0;
     fHistLeadingJetPt[i] = 0;
     fHist2LeadingJetPt[i] = 0;
+    fHistLeadingJetCorrPt[i] = 0;
+    fHistJetPhiEta[i] = 0;
+    fHistJetsPtArea[i] = 0;
+    fHistJetsCorrPtArea[i] = 0;
     fHistJetsNEFvsPt[i] = 0;
     fHistJetsZvsPt[i] = 0;
-    fHistMaxTrackPtvsJetPt[i] = 0;
-    fHistMaxClusPtvsJetPt[i] = 0;
-    fHistMaxPartPtvsJetPt[i] = 0;
-    fHistMaxTrackPtvsJetCorrPt[i] = 0;
-    fHistMaxClusPtvsJetCorrPt[i] = 0;
-    fHistMaxPartPtvsJetCorrPt[i] = 0;
     fHistConstituents[i] = 0;
-    fHistRho[i] = 0;
-    fHistJetsCorrPt[i] = 0;
-    fHistJetsCorrPtArea[i] = 0;
-    fHistLeadingJetCorrPt[i] = 0;
-    fHistRCPtRigid[i] = 0;
-    fHistRCPt[i] = 0;
-    fHistRCPtExLJ[i] = 0;
-    fHistRCPtRand[i] = 0;   
-    fHistDeltaPtRCRigid[i] = 0;
-    fHistDeltaPtRC[i] = 0;
-    fHistDeltaPtRCExLJ[i] = 0;
-    fHistDeltaPtRCRand[i] = 0;
-    fHistEmbNotFoundPhiEta[i] = 0;
-    fHistEmbJetsPt[i] = 0;
-    fHistEmbJetsCorrPt[i] = 0;
-    fHistEmbPartPt[i] = 0;
-    fHistDistEmbPartJetAxis[i] = 0;
-    fHistDeltaPtEmb[i] = 0;
+    fHistTracksJetPt[i] = 0;
+    fHistClustersJetPt[i] = 0;
   }
+
+  SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
@@ -179,63 +88,41 @@ AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
   // Destructor.
 }
 
+//________________________________________________________________________
+Float_t* AliAnalysisTaskSAJF::GenerateFixedBinArray(Int_t n, Float_t min, Float_t max) const
+{
+  Float_t *bins = new Float_t[n+1];
+
+  Float_t binWidth = (max-min)/n;
+  bins[0] = min;
+  for (Int_t i = 1; i <= n; i++) {
+    bins[i] = bins[i-1]+binWidth;
+  }
+
+  return bins;
+}
+
 //________________________________________________________________________
 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
 {
   // Create user output.
 
   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
-  
-  OpenFile(1);
-  fOutput = new TList();
-  fOutput->SetOwner(); 
-  
-  fHistCentrality = new TH1F("fHistCentrality","fHistCentrality", fNbins, 0, 100);
-  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
-  fHistCentrality->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistCentrality);
-
-  fHistDeltaVectorPt = new TH1F("fHistDeltaVectorPt", "fHistDeltaVectorPt", fNbins, -50, 50);
-  fHistDeltaVectorPt->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
-  fHistDeltaVectorPt->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistDeltaVectorPt);
 
   fHistRhoVSleadJetPt = new TH2F("fHistRhoVSleadJetPt","fHistRhoVSleadJetPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
   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", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
-  fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistRCPhiEta);
-
-  fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 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, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
-  fHistRhoVSRCPt->GetXaxis()->SetTitle("A#rho [GeV/c]");
-  fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
-  fOutput->Add(fHistRhoVSRCPt);
-
-  if (!fEmbJetsName.IsNull()) {
-    fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
-    fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
-    fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
-    fOutput->Add(fHistEmbJetPhiEta);
-    
-    fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
-    fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
-    fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
-    fOutput->Add(fHistEmbPartPhiEta);
-    
-    fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
-    fHistRhoVSEmbBkg->GetXaxis()->SetTitle("A#rho [GeV/c]");
-    fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
-    fOutput->Add(fHistRhoVSEmbBkg);
-  }
+  const Int_t nbinsZ = 12;
+  Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
+
+  Float_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
+  Float_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
+  Float_t *binsArea     = GenerateFixedBinArray(30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
+  Float_t *binsEta      = GenerateFixedBinArray(50,-1, 1);
+  Float_t *binsPhi      = GenerateFixedBinArray(101, 0, TMath::Pi() * 2.02);
+  Float_t *bins120      = GenerateFixedBinArray(120, 0, 1.2);
 
   TString histname;
 
@@ -252,283 +139,118 @@ void AliAnalysisTaskSAJF::UserCreateOutputObjects()
     fHistEvents[i]->GetXaxis()->SetBinLabel(5, "OK");
     fOutput->Add(fHistEvents[i]);
 
-    if (fAnaType != kEMCALOnly) { 
-      histname = "fHistTracksPt_";
-      histname += i;
-      fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins / 2.5), fMinBinPt, fMaxBinPt / 2.5);
-      fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-      fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistTracksPt[i]);
-    }
-
-    if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {    
-      histname = "fHistClustersPt_";
-      histname += i;
-      fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(),  (Int_t)(fNbins / 2.5), fMinBinPt, fMaxBinPt / 2.5);
-      fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-      fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistClustersPt[i]);
-    }
-
-    histname = "fHistJetPhiEta_";
-    histname += i;
-    fHistJetPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
-    fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
-    fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
-    fOutput->Add(fHistJetPhiEta[i]);
-
-    histname = "fHistJetsPt_";
-    histname += i;
-    fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
-    fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
-    fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistJetsPt[i]);
-
-    histname = "fHistJetsPtArea_";
-    histname += i;
-    fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
-    fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
-    fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
-    fOutput->Add(fHistJetsPtArea[i]);
-
     histname = "fHistLeadingJetPt_";
     histname += i;
     fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
-    fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV]");
+    fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
     fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistLeadingJetPt[i]);
 
     histname = "fHist2LeadingJetPt_";
     histname += i;
     fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
-    fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV]");
+    fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
     fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHist2LeadingJetPt[i]);
 
+    histname = "fHistLeadingJetCorrPt_";
+    histname += i;
+    fHistLeadingJetCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+    fHistLeadingJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
+    fHistLeadingJetCorrPt[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistLeadingJetCorrPt[i]);
+    
+    histname = "fHistJetPhiEta_";
+    histname += i;
+    fHistJetPhiEta[i] = new TH3F(histname.Data(), histname.Data(), 
+                                50, binsEta, 
+                                101, binsPhi, 
+                                nbinsZ, binsZ);
+    fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
+    fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
+    fHistJetPhiEta[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
+    fOutput->Add(fHistJetPhiEta[i]);
+    
+    histname = "fHistJetsPtArea_";
+    histname += i;
+    fHistJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
+                                 fNbins, binsPt, 
+                                 30, binsArea,
+                                 nbinsZ, binsZ);
+    fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
+    fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
+    fHistJetsPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
+    fOutput->Add(fHistJetsPtArea[i]);
+
     histname = "fHistJetsZvsPt_";
     histname += i;
-    fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+    fHistJetsZvsPt[i] = new TH3F(histname.Data(), histname.Data(), 
+                                120, bins120, 
+                                fNbins, binsPt,
+                                nbinsZ, binsZ);
     fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
-    fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
+    fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
+    fHistJetsZvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
     fOutput->Add(fHistJetsZvsPt[i]);
 
-    if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
+    if (!fCaloName.IsNull()) {
       histname = "fHistJetsNEFvsPt_";
       histname += i;
-      fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+      fHistJetsNEFvsPt[i] = new TH3F(histname.Data(), histname.Data(), 
+                                    120, bins120, 
+                                    fNbins, binsPt,
+                                    nbinsZ, binsZ);
       fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
-      fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
+      fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
+      fHistJetsNEFvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
       fOutput->Add(fHistJetsNEFvsPt[i]);
     }
 
-    if (fAnaType != kEMCALOnly) { 
-      histname = "fHistMaxTrackPtvsJetPt_";
-      histname += i;
-      fHistMaxTrackPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt,  (Int_t)(fNbins / 2.5), fMinBinPt, fMaxBinPt / 2.5);
-      fHistMaxTrackPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
-      fHistMaxTrackPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{track} [GeV/c]");
-      fOutput->Add(fHistMaxTrackPtvsJetPt[i]);
-    }
-
-    if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-      histname = "fHistMaxClusPtvsJetPt_";
-      histname += i;
-      fHistMaxClusPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt,  (Int_t)(fNbins / 2.5), fMinBinPt, fMaxBinPt / 2.5);
-      fHistMaxClusPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
-      fHistMaxClusPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{clus} [GeV/c]");
-      fOutput->Add(fHistMaxClusPtvsJetPt[i]);
-    }
-
-    histname = "fHistMaxPartPtvsJetPt_";
-    histname += i;
-    fHistMaxPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt,  (Int_t)(fNbins / 2.5), fMinBinPt, fMaxBinPt / 2.5);
-    fHistMaxPartPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
-    fHistMaxPartPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{part} [GeV/c]");
-    fOutput->Add(fHistMaxPartPtvsJetPt[i]);
-
-    if (fAnaType != kEMCALOnly) { 
-      histname = "fHistMaxTrackPtvsJetCorrPt_";
-      histname += i;
-      fHistMaxTrackPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt,  (Int_t)(fNbins / 2.5), fMinBinPt, fMaxBinPt / 2.5);
-      fHistMaxTrackPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
-      fHistMaxTrackPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{track} [GeV/c]");
-      fOutput->Add(fHistMaxTrackPtvsJetCorrPt[i]);
-    }
-
-    if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-      histname = "fHistMaxClusPtvsJetCorrPt_";
-      histname += i;
-      fHistMaxClusPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt,  (Int_t)(fNbins / 2.5), fMinBinPt, fMaxBinPt / 2.5);
-      fHistMaxClusPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
-      fHistMaxClusPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{clus} [GeV/c]");
-      fOutput->Add(fHistMaxClusPtvsJetCorrPt[i]);
-    }
-
-    histname = "fHistMaxPartPtvsJetCorrPt_";
-    histname += i;
-    fHistMaxPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, (Int_t)(fNbins / 2.5), fMinBinPt, fMaxBinPt / 2.5);
-    fHistMaxPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
-    fHistMaxPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{part} [GeV/c]");
-    fOutput->Add(fHistMaxPartPtvsJetCorrPt[i]);
-
-    histname = "fHistConstituents_";
-    histname += i;
-    fHistConstituents[i] = new TH2F(histname.Data(), histname.Data(), 100, 1, 101, 100, -0.5, 99.5);
-    fHistConstituents[i]->GetXaxis()->SetTitle("p_{T,part} [GeV/c]");
-    fHistConstituents[i]->GetYaxis()->SetTitle("no. of particles");
-    fOutput->Add(fHistConstituents[i]);
-
-    histname = "fHistRho_";
-    histname += i;
-    fHistRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
-    fHistRho[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-    fHistRho[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRho[i]);
-
-    histname = "fHistJetsCorrPt_";
-    histname += i;
-    fHistJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistJetsCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
-    fHistJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistJetsCorrPt[i]);
-
     histname = "fHistJetsCorrPtArea_";
     histname += i;
-    fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
+    fHistJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
+                                     fNbins * 2, binsCorrPt, 
+                                     30, binsArea,
+                                     nbinsZ, binsZ);
     fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
     fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
     fOutput->Add(fHistJetsCorrPtArea[i]);
 
-    histname = "fHistLeadingJetCorrPt_";
-    histname += i;
-    fHistLeadingJetCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistLeadingJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{RC} [GeV/c]");
-    fHistLeadingJetCorrPt[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistLeadingJetCorrPt[i]);
-
-    histname = "fHistRCPtRigid_";
-    histname += i;
-    fHistRCPtRigid[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
-    fHistRCPtRigid[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
-    fHistRCPtRigid[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRCPtRigid[i]);
-
-    histname = "fHistRCPt_";
-    histname += i;
-    fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 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, fMinBinPt, fMaxBinPt * 2);
-    fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
-    fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRCPtExLJ[i]);
-
-    histname = "fHistRCPtRand_";
-    histname += i;
-    fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
-    fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
-    fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRCPtRand[i]);
-
-    histname = "fHistDeltaPtRCRigid_";
-    histname += i;
-    fHistDeltaPtRCRigid[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistDeltaPtRCRigid[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
-    fHistDeltaPtRCRigid[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPtRCRigid[i]);
-
-    histname = "fHistDeltaPtRC_";
-    histname += i;
-    fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
-    fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPtRC[i]);
-
-    histname = "fHistDeltaPtRCExLJ_";
+    histname = "fHistConstituents_";
     histname += i;
-    fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
-    fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPtRCExLJ[i]);
+    fHistConstituents[i] = new TH3F(histname.Data(), histname.Data(), 100, 1, 101, 100, -0.5, 99.5, fNbins, fMinBinPt, fMaxBinPt);
+    fHistConstituents[i]->GetXaxis()->SetTitle("p_{T,part} (GeV/c)");
+    fHistConstituents[i]->GetYaxis()->SetTitle("no. of particles");
+    fHistConstituents[i]->GetZaxis()->SetTitle("p_{T,jet} (GeV/c)");
+    fOutput->Add(fHistConstituents[i]);
 
-    histname = "fHistDeltaPtRCRand_";
+    histname = "fHistTracksJetPt_";
     histname += i;
-    fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
-    fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPtRCRand[i]);
-
-    if (!fEmbJetsName.IsNull()) {
-      histname = "fHistEmbJetsPt_";
+    fHistTracksJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
+    fHistTracksJetPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
+    fHistTracksJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
+    fHistTracksJetPt[i]->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistTracksJetPt[i]);
+
+    if (!fCaloName.IsNull()) {
+      histname = "fHistClustersJetPt_";
       histname += i;
-      fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
-      fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
-      fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistEmbJetsPt[i]);
-
-      histname = "fHistEmbJetsCorrPt_";
-      histname += i;
-      fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-      fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
-      fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistEmbJetsCorrPt[i]);
-
-      histname = "fHistEmbJetsArea_";
-      histname += i;
-      fHistEmbJetsArea[i] = new TH1F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
-      fHistEmbJetsArea[i]->GetXaxis()->SetTitle("area");
-      fHistEmbJetsArea[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistEmbJetsArea[i]);
-
-      histname = "fHistEmbPartPt_";
-      histname += i;
-      fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
-      fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
-      fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistEmbPartPt[i]);
-
-      histname = "fHistDistEmbPartJetAxis_";
-      histname += i;
-      fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
-      fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
-      fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistDistEmbPartJetAxis[i]);
-
-      histname = "fHistEmbNotFoundPhiEta_";
-      histname += i;
-      fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
-      fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
-      fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
-      fOutput->Add(fHistEmbNotFoundPhiEta[i]);
-      
-      histname = "fHistDeltaPtEmb_";
-      histname += i;
-      fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-      fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [GeV/c]");
-      fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistDeltaPtEmb[i]);
+      fHistClustersJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
+      fHistClustersJetPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
+      fHistClustersJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
+      fHistClustersJetPt[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistClustersJetPt[i]);
     }
   }
 
   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
-}
 
-//________________________________________________________________________
-Bool_t AliAnalysisTaskSAJF::RetrieveEventObjects()
-{
-  // Retrieve event objects.
-
-  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
-    return kFALSE;
-  
-  if (fRho)
-    fRhoVal = fRho->GetVal();
-
-  return kTRUE;
+  delete binsPt;
+  delete binsCorrPt;
+  delete binsArea;
+  delete binsEta;
+  delete binsPhi;
+  delete bins120;
 }
 
 //________________________________________________________________________
@@ -536,278 +258,60 @@ Bool_t AliAnalysisTaskSAJF::FillHistograms()
 {
   // Fill histograms.
 
-  Int_t maxJetIndex  = -1;
-  Int_t max2JetIndex = -1;
+  if (!fJets) {
+    AliError(Form("%s - Jet array not provided, returning...", GetName()));
+    return kFALSE;
+  }
+
+  if (fJets->GetEntriesFast() < 1) { // no jets in array, skipping
+    fHistEvents[fCentBin]->Fill("No jets", 1);
+    return kTRUE;
+  }
 
-  GetLeadingJets(maxJetIndex, max2JetIndex);
+  Int_t *sortedJets = GetSortedArray(fJets);
   
-  if (maxJetIndex < 0) { // no accepted jet, skipping
+  if (!sortedJets || sortedJets[0] < 0) { // no accepted jets, skipping
     fHistEvents[fCentBin]->Fill("No jets", 1);
     return kTRUE;
   }
 
-  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(maxJetIndex));
-  if (!jet) {  // error, I cannot get the lead jet from collection (should never happen), skipping
+  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
+  if (!jet) {  // error, I cannot get the leading jet from collection (should never happen), skipping
     fHistEvents[fCentBin]->Fill("Max jet not found", 1);
     return kTRUE;
   }
 
   // OK, event accepted
 
+  Float_t maxJetCorrPt = jet->Pt() - fRhoVal * jet->Area();
+
   if (fRhoVal == 0) 
     fHistEvents[fCentBin]->Fill("Rho == 0", 1);
 
-  Float_t maxJetCorrPt = jet->Pt() - fRhoVal * jet->Area();
-
-  if (maxJetCorrPt <= 0)
+  else if (maxJetCorrPt <= 0)
     fHistEvents[fCentBin]->Fill("Max jet <= 0", 1);
 
-  fHistEvents[fCentBin]->Fill("OK", 1);
-
-  // ************
-  // General histograms
-  // _________________________________
-
-  DoJetLoop();
-  DoTrackLoop();
-  DoClusterLoop();
-
-  fHistCentrality->Fill(fCent);
-  fHistRho[fCentBin]->Fill(fRhoVal);
+  else
+    fHistEvents[fCentBin]->Fill("OK", 1);
 
   if (jet) {
     fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
     fHistRhoVSleadJetPt->Fill(fRhoVal * jet->Area(), jet->Pt());
     fHistLeadingJetCorrPt[fCentBin]->Fill(maxJetCorrPt);
   }
-  
+
   AliEmcalJet* jet2 = 0;
-  if (max2JetIndex >= 0)
-    jet2 = static_cast<AliEmcalJet*>(fJets->At(max2JetIndex));
+  if (sortedJets[1] >= 0)
+    jet2 = static_cast<AliEmcalJet*>(fJets->At(sortedJets[1]));
 
   if (jet2)
     fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
 
-  // ************
-  // Random cones
-  // _________________________________
-  
-  const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
-
-  // Simple random cones
-  for (Int_t i = 0; i < fRCperEvent; i++) {
-    Float_t RCpt = 0;
-    Float_t RCptRigid = 0;
-    Float_t RCeta = 0;
-    Float_t RCphi = 0;
-    GetRigidCone(RCpt, RCptRigid, RCeta, RCphi, 0);
-    if (RCpt > 0) {
-      fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
-      fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
-    }
-    if (RCptRigid > 0) {
-      fHistRCPtRigid[fCentBin]->Fill(RCptRigid / rcArea);
-      fHistDeltaPtRCRigid[fCentBin]->Fill(RCptRigid - rcArea * fRhoVal);
-    }
-  
-    // Random cones far from leading jet
-    RCpt = 0;
-    RCptRigid = 0;
-    RCeta = 0;
-    RCphi = 0;
-    GetRigidCone(RCpt, RCptRigid, RCeta, RCphi, jet);
-    if (RCpt > 0) {
-      fHistRCPhiEta->Fill(RCeta, RCphi);
-      fHistRhoVSRCPt->Fill(fRhoVal, RCpt / rcArea);
-      
-      Float_t dphi = RCphi - jet->Phi();
-      if (dphi > 4.8) dphi -= TMath::Pi() * 2;
-      if (dphi < -1.6) dphi += TMath::Pi() * 2; 
-      fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
-    
-      fHistRCPtExLJ[fCentBin]->Fill(RCpt / rcArea);
-      fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
-    }
-    
-    // Random cones with randomized particles
-    RCpt = 0;
-    RCptRigid = 0;
-    RCeta = 0;
-    RCphi = 0;
-    GetRigidCone(RCpt, RCptRigid, RCeta, RCphi, 0, fRandTracks, fRandCaloClusters);
-    if (RCpt > 0) {
-      fHistRCPtRand[fCentBin]->Fill(RCpt / rcArea);
-      fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
-    }  
-  }
-  
-  // ************
-  // Embedding
-  // _________________________________
-
-  if (!fEmbJets)
-    return kTRUE;
-
-  AliEmcalJet *embJet  = 0;
-  TObject     *embPart = 0;
-
-  DoEmbJetLoop(embJet, embPart);
-
-  if (embJet) {
-    Double_t probePt = 0;
-    Double_t probeEta = 0;
-    Double_t probePhi = 0;
-
-    AliVCluster *cluster = dynamic_cast<AliVCluster*>(embPart);
-    if (cluster) {
-      TLorentzVector nPart;
-      cluster->GetMomentum(nPart, fVertex);
-
-      probeEta = nPart.Eta();
-      probePhi = nPart.Phi();
-      probePt = nPart.Pt();
-    }
-    else {
-      AliVTrack *track = dynamic_cast<AliVTrack*>(embPart);
-      if (track) {
-       probeEta = track->Eta();
-       probePhi = track->Phi();
-       probePt = track->Pt();
-      }
-      else {
-       AliWarning(Form("%s - Embedded jet found but embedded particle not found (wrong type?)!", GetName()));
-       return kTRUE;
-      }
-    }
-    Double_t distProbeJet = TMath::Sqrt((embJet->Eta() - probeEta) * (embJet->Eta() - probeEta) + (embJet->Phi() - probePhi) * (embJet->Phi() - probePhi));
-
-    fHistEmbPartPt[fCentBin]->Fill(probePt);
-    fHistEmbPartPhiEta->Fill(probeEta, probePhi);
-    
-    fHistDistEmbPartJetAxis[fCentBin]->Fill(distProbeJet);
-
-    fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
-    fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - fRhoVal * embJet->Area());
-    fHistEmbJetsArea[fCentBin]->Fill(embJet->Area());
-    fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
-
-    fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRhoVal - probePt);
-    fHistRhoVSEmbBkg->Fill(embJet->Area() * fRhoVal, embJet->Pt() - probePt);
-  }
-  else {
-    if (fAnaType != kEMCALOnly)
-      DoEmbTrackLoop();
-    if (fAnaType == kEMCAL || fAnaType == kEMCALOnly)
-      DoEmbClusterLoop();
-    if (fEmbeddedTrackId >= 0) {
-      AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
-      fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
-    }
-    else if (fEmbeddedClusterId >= 0) {
-      AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
-      TLorentzVector nPart;
-      cluster2->GetMomentum(nPart, fVertex);
-      fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
-    }
-    else {
-      AliWarning(Form("%s - Embedded particle not found!", GetName()));
-    }
-  }
+  DoJetLoop();
 
   return kTRUE;
 }
 
-//________________________________________________________________________
-void AliAnalysisTaskSAJF::GetLeadingJets(Int_t &maxJetIndex, Int_t &max2JetIndex)
-{
-  // Get the leading jets.
-
-  if (!fJets)
-    return;
-
-  const Int_t njets = fJets->GetEntriesFast();
-
-  Float_t maxJetPt = -999;
-  Float_t max2JetPt = -999;
-  for (Int_t ij = 0; ij < njets; ij++) {
-
-    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
-
-    if (!jet) {
-      AliError(Form("Could not receive jet %d", ij));
-      continue;
-    }  
-
-    if (!AcceptJet(jet))
-      continue;
-
-    Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
-
-    if (maxJetIndex == -1 || maxJetPt < corrPt) {
-      max2JetPt = maxJetPt;
-      max2JetIndex = maxJetIndex;
-      maxJetPt = corrPt;
-      maxJetIndex = ij;
-    }
-    else if (max2JetIndex == -1 || max2JetPt < corrPt) {
-      max2JetPt = corrPt;
-      max2JetIndex = ij;
-    }
-  }
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSAJF::DoClusterLoop()
-{
-  // Do cluster loop.
-
-  if (!fCaloClusters)
-    return;
-
-  Int_t nclusters =  fCaloClusters->GetEntriesFast();
-
-  for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
-    AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
-    if (!cluster) {
-      AliError(Form("Could not receive cluster %d", iClusters));
-      continue;
-    }  
-
-    if (!AcceptCluster(cluster)) 
-      continue;
-
-    fHistClustersPt[fCentBin]->Fill(cluster->E());
-  }
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSAJF::DoTrackLoop()
-{
-  // Do track loop.
-
-  if (!fTracks)
-    return;
-
-  Int_t ntracks = fTracks->GetEntriesFast();
-
-  for (Int_t i = 0; i < ntracks; i++) {
-
-    AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track  
-
-    if (!track) {
-      AliError(Form("Could not retrieve track %d",i)); 
-      continue; 
-    }
-
-    AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
-    
-    if (vtrack && !AcceptTrack(vtrack)) 
-      continue;
-    
-    fHistTracksPt[fCentBin]->Fill(track->Pt());
-  }
-}
-
 //________________________________________________________________________
 void AliAnalysisTaskSAJF::DoJetLoop()
 {
@@ -835,36 +339,29 @@ void AliAnalysisTaskSAJF::DoJetLoop()
 
     Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
 
-    fHistJetsPt[fCentBin]->Fill(jet->Pt());
-    fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
-    fHistJetsCorrPt[fCentBin]->Fill(corrPt);
-    fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
-
-    fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
+    Float_t ptLeading = 0;
 
-    fHistMaxPartPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxPartPt());
-    fHistMaxPartPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxPartPt());
+    if (fLeadingHadronType == 0)       // charged leading hadron
+      ptLeading = jet->MaxTrackPt();
+    else if (fLeadingHadronType == 1)  // neutral leading hadron
+      ptLeading = jet->MaxClusterPt();
+    else                               // charged or neutral
+      ptLeading = jet->MaxPartPt();
 
-    if (fAnaType != kEMCALOnly) {
-      fHistMaxTrackPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxTrackPt());
-      fHistMaxTrackPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxTrackPt());
-    }
-
-    if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-      fHistMaxClusPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxClusterPt());
-      fHistMaxClusPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxClusterPt());
-      fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt());
-    }
+    fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi(), ptLeading);
+    fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area(), ptLeading);
+    fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area(), ptLeading);
 
-    Float_t scalarpt = 0;
+    if (fCaloClusters)
+      fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt(), ptLeading);
 
     if (fTracks) {
       for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
        AliVParticle *track = jet->TrackAt(it, fTracks);
        if (track) {
-         fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt());
+         fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt(), ptLeading);
          constituents.Fill(track->Pt());
-         scalarpt += track->Pt();
+         fHistTracksJetPt[fCentBin]->Fill(track->Pt(), jet->Pt());
        }
       }
     }
@@ -876,372 +373,21 @@ void AliAnalysisTaskSAJF::DoJetLoop()
        if (cluster) {
          TLorentzVector nPart;
          cluster->GetMomentum(nPart, fVertex);
-         fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt());
-         scalarpt += nPart.Pt();
+         fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt(), ptLeading);
          constituents.Fill(nPart.Pt());
+         fHistClustersJetPt[fCentBin]->Fill(nPart.Pt(), jet->Pt());
        }
       }
     }
 
-    fHistDeltaVectorPt->Fill(scalarpt - jet->Pt());
-
     for (Int_t i = 1; i <= constituents.GetNbinsX(); i++) {
-      fHistConstituents[fCentBin]->Fill(constituents.GetBinCenter(i), constituents.GetBinContent(i));
+      fHistConstituents[fCentBin]->Fill(constituents.GetBinCenter(i), constituents.GetBinContent(i), jet->Pt());
     }
 
     constituents.Reset();
   } //jet loop 
 }
 
-//________________________________________________________________________
-void AliAnalysisTaskSAJF::DoEmbTrackLoop()
-{
-  // Do track loop.
-
-  if (!fEmbTracks)
-    return;
-
-  fEmbeddedTrackId = -1;
-
-  Int_t ntracks = fEmbTracks->GetEntriesFast();
-
-  for (Int_t i = 0; i < ntracks; i++) {
-
-    AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track  
-
-    if (!track) {
-      AliError(Form("Could not retrieve track %d",i)); 
-      continue; 
-    }
-
-    AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
-    
-    if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
-      continue;
-
-    if (track->GetLabel() == 100) {
-      fEmbeddedTrackId = i;
-      continue;
-    }
-  }
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSAJF::DoEmbClusterLoop()
-{
-  // Do cluster loop.
-
-  if (!fEmbCaloClusters)
-    return;
-
-  fEmbeddedClusterId = -1;
-
-  Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
-
-  for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
-    AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
-    if (!cluster) {
-      AliError(Form("Could not receive cluster %d", iClusters));
-      continue;
-    }  
-
-    if (!AcceptCluster(cluster, kTRUE)) 
-      continue;
-
-    if (cluster->Chi2() == 100) {
-      fEmbeddedClusterId = iClusters;
-      continue;
-    }
-  }
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
-{
-  // Do the embedded jet loop.
-
-  if (!fEmbJets)
-    return;
-
-  TLorentzVector maxClusVect;
-
-  Int_t nembjets = fEmbJets->GetEntriesFast();
-
-  for (Int_t ij = 0; ij < nembjets; ij++) {
-      
-    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
-      
-    if (!jet) {
-      AliError(Form("Could not receive jet %d", ij));
-      continue;
-    } 
-      
-    if (!AcceptJet(jet))
-      continue;
-
-    if (!jet->IsMC())
-      continue;
-
-    AliVParticle *maxTrack = 0;
-
-    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
-      AliVParticle *track = jet->TrackAt(it, fEmbTracks);
-      
-      if (!track) 
-       continue;
-
-      if (track->GetLabel() != 100)
-       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, fEmbCaloClusters);
-      
-      if (!cluster) 
-       continue;
-
-      if (cluster->Chi2() != 100)
-       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)) {
-      embPart = maxClus;
-      embJet = jet;
-    }
-    else if (maxTrack) {
-      embPart = maxTrack;
-      embJet = jet;
-    }
-
-    return;  // MC jet found, exit
-  }
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSAJF::GetRigidCone(Float_t &pt, Float_t &ptrigid, Float_t &eta, Float_t &phi,
-                                      AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
-{
-  // Get rigid cone.
-
-  if (!tracks)
-    tracks = fTracks;
-
-  if (!clusters)
-    clusters = fCaloClusters;
-
-  eta = 0;
-  phi = 0;
-  pt = 0;
-  ptrigid = 0;
-
-  if (!tracks && !clusters)
-    return;
-
-  Float_t LJeta = 999;
-  Float_t LJphi = 999;
-
-  if (jet) {
-    LJeta = jet->Eta();
-    LJphi = jet->Phi();
-  }
-
-  Float_t maxEta = fMaxEta;
-  Float_t minEta = fMinEta;
-  Float_t maxPhi = fMaxPhi;
-  Float_t minPhi = fMinPhi;
-
-  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) {
-    AliWarning(Form("%s: Could not get random cone!", GetName()));
-    return;
-  }
-
-  TVector3 rigidAxis;
-  rigidAxis.SetPtEtaPhi(1, eta, phi);
-  rigidAxis = rigidAxis.Unit();
-
-  if ((fAnaType == kEMCAL || fAnaType == kEMCALOnly) && clusters) {
-    Int_t nclusters =  clusters->GetEntriesFast();
-    for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
-      AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
-      if (!cluster) {
-       AliError(Form("Could not receive cluster %d", iClusters));
-       continue;
-      }  
-      
-      if (!AcceptCluster(cluster, fMCAna))
-       continue;
-      
-      TLorentzVector nPart;
-      cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
-     
-      Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
-
-      if (d <= fJetRadius) {
-       TVector3 vect = nPart.Vect();
-       vect *= vect * rigidAxis / vect.Mag();
-       ptrigid += vect.Pt();
-       
-       pt += nPart.Pt();
-      }
-    }
-  }
-
-  if (tracks) {
-    Int_t ntracks = tracks->GetEntriesFast();
-    for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
-      AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));         
-      if(!track) {
-       AliError(Form("Could not retrieve track %d",iTracks)); 
-       continue; 
-      }
-      
-      if (!AcceptTrack(track, fMCAna)) 
-       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){
-       TVector3 vect(track->Px(), track->Py(), track->Pz());
-       vect *= vect * rigidAxis / vect.Mag();
-       ptrigid += vect.Pt();
-       
-       pt += track->Pt();
-      }
-    }
-  }
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSAJF::ExecOnce()
-{
-  // Initialize the analysis.
-
-  if (!fRhoName.IsNull() && !fRho) {
-    fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
-    if (!fRho) {
-      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
-      return;
-    }
-  }
-  
-  if (!fEmbJetsName.IsNull() && !fEmbJets) {
-    fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
-    if (!fEmbJets) {
-      AliError(Form("%s: Could not retrieve emb jets %s!", GetName(), fEmbJetsName.Data()));
-      return;
-    }
-    else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
-      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data())); 
-      fEmbJets = 0;
-      return;
-    }
-  }
-
-  if (!fEmbCaloName.IsNull() && (fAnaType == kEMCAL || fAnaType == kEMCALOnly) && !fEmbCaloClusters) {
-    fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
-    if (!fEmbCaloClusters) {
-      AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
-      return;
-    }
-    else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
-      AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data())); 
-      fEmbCaloClusters = 0;
-      return;
-    }
-  }
-
-  if (!fEmbTracksName.IsNull() && fAnaType != kEMCALOnly && !fEmbTracks) {
-    fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
-    if (!fEmbTracks) {
-      AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
-      return;
-    }
-    else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
-      AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data())); 
-      fEmbTracks = 0;
-      return;
-    }
-  }
-
-  if (!fRandCaloName.IsNull() && (fAnaType == kEMCAL || fAnaType == kEMCALOnly) && !fRandCaloClusters) {
-    fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
-    if (!fRandCaloClusters) {
-      AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
-      return;
-    }
-    else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
-      AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data())); 
-      fRandCaloClusters = 0;
-      return;
-    }
-  }
-
-  if (!fRandTracksName.IsNull() && fAnaType != kEMCALOnly && !fRandTracks) {
-    fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
-    if (!fRandTracks) {
-      AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
-      return;
-    }
-    else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
-      AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data())); 
-      fRandTracks = 0;
-      return;
-    }
-  }
-
-  AliAnalysisTaskEmcalJet::ExecOnce();
-
-  if (fRCperEvent < 0) {
-    Double_t area = (fMaxEta - fMinEta) * (fMaxPhi - fMinPhi);
-    Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
-    fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
-    if (fRCperEvent == 0)
-      fRCperEvent = 1;
-  }
-
-  if (fMinRC2LJ < 0)
-    fMinRC2LJ = fJetRadius * 1.5;
-
-  const Float_t maxDist = TMath::Max(fMaxPhi - fMinPhi, fMaxEta - fMinEta) / 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;
-  }
-}
-
 //________________________________________________________________________
 void AliAnalysisTaskSAJF::Terminate(Option_t *) 
 {
index 18d6a2e4c5d70a29b27a74272b8c797c5aa5b0bf..7ffba44e2f407d077b17cfdc372e210b9e1c7ccc 100644 (file)
@@ -7,6 +7,7 @@ class TClonesArray;
 class TString;
 class TH1F;
 class TH2F;
+class TH3F;
 class AliRhoParameter;
 
 #include "AliAnalysisTaskEmcalJet.h"
@@ -21,107 +22,36 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskEmcalJet {
   void                        UserCreateOutputObjects();
   void                        Terminate(Option_t *option);
 
-  void                        SetJetMinRC2LJ(Float_t d)                            { fMinRC2LJ                = d          ; } 
-  void                        SetEmbJetsName(const char *n)                        { fEmbJetsName             = n          ; }
-  void                        SetRandTracksName(const char *n)                     { fRandTracksName          = n          ; }
-  void                        SetRandClusName(const char *n)                       { fRandCaloName            = n          ; }
-  void                        SetEmbTracksName(const char *n)                      { fEmbTracksName           = n          ; }
-  void                        SetEmbClusName(const char *n)                        { fEmbCaloName             = n          ; }
-  void                        SetRhoName(const char *n)                            { fRhoName                 = n          ; }
-  void                        SetMC(Bool_t m)                                      { fMCAna                   = m          ; }
-  void                        SetRCperEvent(Int_t n)                               { fRCperEvent              = n          ; }
+  void                        SetLeadingHadronType(Int_t t)     { fLeadingHadronType = t; }
 
  protected:
-  void                        ExecOnce()                                                                                    ;
-  Bool_t                      RetrieveEventObjects()                                                                        ;
-  Bool_t                      FillHistograms()                                                                              ;
-  void                        GetLeadingJets(Int_t &maxJetIndex, Int_t &max2JetIndex)                                       ;
-  void                        DoJetLoop()                                                                                   ;
-  void                        DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)                                         ;
-  void                        DoTrackLoop()                                                                                 ;
-  void                        DoClusterLoop()                                                                               ;
-  void                        DoEmbTrackLoop()                                                                              ;
-  void                        DoEmbClusterLoop()                                                                            ;
-  void                        GetRigidCone(Float_t &pt, Float_t &ptrigid, Float_t &eta, Float_t &phi, 
-                                          AliEmcalJet *jet = 0, TClonesArray* tracks = 0, TClonesArray* clusters = 0) const;
+  Bool_t                      FillHistograms()                                              ;
+  void                        DoJetLoop()                                                   ;
+  Float_t*                    GenerateFixedBinArray(Int_t n, Float_t min, Float_t max) const;
 
-  Bool_t                      fMCAna;                      // =true MC analysis (toy model)
-  Float_t                     fMinRC2LJ;                   // Minimum distance random cone to leading jet
-  TString                     fEmbJetsName;                // Name of embedded jet collection
-  TString                     fEmbTracksName;              // Name of embedded track collection
-  TString                     fEmbCaloName;                // Name of embedded calo cluster collection
-  TString                     fRandTracksName;             // Name of randomized track collection
-  TString                     fRandCaloName;               // Name of randomized calo cluster collection
-  TString                     fRhoName;                    // Name of rho object
-  Int_t                       fRCperEvent;                 // No. of random cones per event
-
-  TClonesArray               *fEmbJets;                    //!Embedded Jets
-  TClonesArray               *fEmbTracks;                  //!Embedded tracks
-  TClonesArray               *fEmbCaloClusters;            //!Embedded clusters  
-  TClonesArray               *fRandTracks;                 //!Randomized tracks
-  TClonesArray               *fRandCaloClusters;           //!Randomized clusters
-  AliRhoParameter            *fRho;                        //!Event rho
-  Double_t                    fRhoVal;                     //!Event rho value
-  Int_t                       fEmbeddedClusterId;          //!Embedded cluster id
-  Int_t                       fEmbeddedTrackId;            //!Embedded track id
+  Int_t                       fLeadingHadronType;          // 0 = charged, 1 = neutral, 2 = both
 
   // General histograms
-  TH1F                       *fHistCentrality;             //!Event centrality distribution
   TH1F                       *fHistEvents[4];              //!Events accepted/rejected
-  TH1F                       *fHistTracksPt[4];            //!Inclusive track pt spectrum
-  TH1F                       *fHistClustersPt[4];          //!Inclusive clusters pt spectrum
-  TH2F                       *fHistJetPhiEta[4];           //!Phi-Eta distribution of jets
-  TH1F                       *fHistJetsPt[4];              //!Inclusive jet pt spectrum
-  TH2F                       *fHistJetsPtArea[4];          //!Jet pt vs. area
   TH1F                       *fHistLeadingJetPt[4];        //!Leading jet pt spectrum
   TH1F                       *fHist2LeadingJetPt[4];       //!Second leading jet pt spectrum
-  TH2F                       *fHistJetsNEFvsPt[4];         //!Jet neutral energy fraction vs. jet pt
-  TH2F                       *fHistJetsZvsPt[4];           //!Constituent Pt over Jet Pt ratio vs. jet pt
-  TH2F                       *fHistMaxTrackPtvsJetPt[4];   //!Max constituent track pt vs. jet pt
-  TH2F                       *fHistMaxClusPtvsJetPt[4];    //!Max constituent cluster pt vs. jet pt
-  TH2F                       *fHistMaxPartPtvsJetPt[4];    //!Max constituent particle (track or cluster) pt vs. jet pt
-  TH2F                       *fHistMaxTrackPtvsJetCorrPt[4];   //!Max constituent track pt vs. jet pt
-  TH2F                       *fHistMaxClusPtvsJetCorrPt[4];    //!Max constituent cluster pt vs. jet pt
-  TH2F                       *fHistMaxPartPtvsJetCorrPt[4];    //!Max constituent particle (track or cluster) pt vs. jet pt
-  TH2F                       *fHistConstituents[4];            //!x axis = constituents pt; y axis = no. of constituents; z axis = jet pt
-  TH1F                       *fHistDeltaVectorPt;              //!Delta Pt between vector and scalar sum
-
-  // Rho
-  TH1F                       *fHistRho[4];                    //!Rho distribution
-  TH2F                       *fHistRhoVSleadJetPt;            //!Area(leadjetarea) * rho vs. leading jet pt
-  TH1F                       *fHistJetsCorrPt[4];             //!Corrected inclusive jet pt spectrum
-  TH2F                       *fHistJetsCorrPtArea[4];         //!Jet pt vs. area
-  TH1F                       *fHistLeadingJetCorrPt[4];       //!Corrected leading jet pt spectrum
-
-  // Random cones
-  TH2F                       *fHistRCPhiEta;               //!Phi-Eta distribution of random cones
-  TH1F                       *fHistRCPtRigid[4];           //!Random cone pt, rigid
-  TH1F                       *fHistRCPt[4];                //!Random cone pt
-  TH1F                       *fHistRCPtExLJ[4];            //!Random cone pt, imposing min distance from leading jet
-  TH1F                       *fHistRCPtRand[4];            //!Random cone pt, randomized particles
-  TH2F                       *fHistRCPtExLJVSDPhiLJ;       //!Random cone pt, imposing min distance from leading jet, vs. deltaPhi leading jet
-  TH2F                       *fHistRhoVSRCPt;              //!Rho vs. Pt(RCExLJ) / Area(RCExLJ)
-  TH1F                       *fHistDeltaPtRCRigid[4];      //!deltaPt = Pt(RC) - A * rho, rigid
-  TH1F                       *fHistDeltaPtRC[4];           //!deltaPt = Pt(RC) - A * rho
-  TH1F                       *fHistDeltaPtRCExLJ[4];       //!deltaPt = Pt(RC) - A * rho, imposing min distance from leading jet
-  TH1F                       *fHistDeltaPtRCRand[4];       //!deltaPt = Pt(RC) - A * rho, randomzied particles
-
-  // Jet embedding
-  TH2F                       *fHistEmbNotFoundPhiEta[4];   //!Phi-Eta of "not found" embedded particles
-  TH1F                       *fHistEmbJetsPt[4];           //!Pt distribution of embedded jets
-  TH1F                       *fHistEmbJetsCorrPt[4];       //!Pt distribution of embedded jets
-  TH1F                       *fHistEmbJetsArea[4];         //!Area distribution of embedded jets
-  TH1F                       *fHistEmbPartPt[4];           //!Pt distribution of embedded particle
-  TH2F                       *fHistEmbJetPhiEta;           //!Phi-Eta distribution of embedded jets
-  TH2F                       *fHistEmbPartPhiEta;          //!Phi-Eta distribution of embedded particles
-  TH1F                       *fHistDistEmbPartJetAxis[4];  //!Distance between embedded particle and jet axis
-  TH2F                       *fHistRhoVSEmbBkg;            //!Area(embjet) * rho vs. Pt(embjet) - Pt(embtrack)
-  TH1F                       *fHistDeltaPtEmb[4];          //!deltaPt = Pt(embjet) - Area(embjet) * rho - Pt(embtrack)
+  TH1F                       *fHistLeadingJetCorrPt[4];    //!Corrected leading jet pt spectrum
+  TH2F                       *fHistRhoVSleadJetPt;         //!Area(leadjet) * rho vs. leading jet pt
+
+  // Inclusive jets histograms
+  TH3F                       *fHistJetPhiEta[4];           //!Phi-Eta distribution of jets
+  TH3F                       *fHistJetsPtArea[4];          //!Jet pt vs. area
+  TH3F                       *fHistJetsCorrPtArea[4];      //!Jet corr pt vs. area
+  TH3F                       *fHistJetsNEFvsPt[4];         //!Jet neutral energy fraction vs. jet pt
+  TH3F                       *fHistJetsZvsPt[4];           //!Constituent Pt over Jet Pt ratio vs. jet pt
+  TH3F                       *fHistConstituents[4];        //!x axis = constituents pt; y axis = no. of constituents; z axis = jet pt
+  TH2F                       *fHistTracksJetPt[4];         //!Track pt vs. jet pt
+  TH2F                       *fHistClustersJetPt[4];       //!Cluster pt vs. jet pt
 
  private:
   AliAnalysisTaskSAJF(const AliAnalysisTaskSAJF&);            // not implemented
   AliAnalysisTaskSAJF &operator=(const AliAnalysisTaskSAJF&); // not implemented
 
-  ClassDef(AliAnalysisTaskSAJF, 10) // jet analysis task
+  ClassDef(AliAnalysisTaskSAJF, 11) // jet analysis task
 };
 #endif
index 158fa4338c283977eb9a91f0e36ff45b5c452c85..fe46d53539a249a765b6302da5240de1cfaf2546 100644 (file)
@@ -32,15 +32,9 @@ ClassImp(AliAnalysisTaskSAQA)
 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() : 
   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
   fCellEnergyCut(0.1),
-  fDoTrigger(kFALSE),
-  fRepropagateTracks(kFALSE),
-  fTrgClusName("ClustersL1GAMMAFEE"),
-  fTrgClusters(0),
   fNclusters(0),
   fNtracks(0),
   fNjets(0),
-  fHistCentrality(0),
-  fHistZVertex(0),
   fHistTracksCent(0),
   fHistClusCent(0),
   fHistJetsCent(0),
@@ -48,18 +42,10 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
   fHistJetsParts(0),
   fHistCellsCent(0),
   fHistCellsTracks(0),
-  fHistMaxL1FastORCent(0),
-  fHistMaxL1ClusCent(0),
-  fHistMaxL1ThrCent(0),
-  fHistTracksPt(0),
-  fHistTrPhiEta(0),
   fHistTrEmcPhiEta(0),
   fHistTrPhiEtaNonProp(0),
   fHistDeltaEtaPt(0),
   fHistDeltaPhiPt(0),
-  fHistDeltaEtaNewProp(0),
-  fHistDeltaPhiNewProp(0),
-  fHistClusPhiEtaEnergy(0),
   fHistNCellsEnergy(0),
   fHistClusTimeEnergy(0),
   fHistCellsEnergy(0),
@@ -74,35 +60,23 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
     fHistTrackEta[i] = 0;
   }
 
-  for (Int_t i = 0; i < 6; i++) {
-    fHistTrackPhiPt[i] = 0;
-    fHistTrackEtaPt[i] = 0;
-  }
-
   for (Int_t i = 0; i < 4; i++) {
-    fHistJetsPhiEta[i] = 0;
-    fHistJetsPtNonBias[i] = 0;
-    fHistJetsPtTrack[i] = 0;
-    fHistJetsPtClus[i] = 0;
-    fHistJetsPt[i] = 0;
+    fHistTrPhiEtaPt[i] = 0;
+    fHistClusPhiEtaEnergy[i] = 0;
+    fHistJetsPhiEtaPt[i] = 0;
     fHistJetsPtArea[i] = 0;
-    fHistJetsPtAreaNonBias[i] = 0;
   }
+
+  SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) : 
   AliAnalysisTaskEmcalJet(name, kTRUE),
   fCellEnergyCut(0.1),
-  fDoTrigger(kFALSE),
-  fRepropagateTracks(kFALSE),
-  fTrgClusName("ClustersL1GAMMAFEE"),
-  fTrgClusters(0),
   fNclusters(0),
   fNtracks(0),
   fNjets(0),
-  fHistCentrality(0),
-  fHistZVertex(0),
   fHistTracksCent(0),
   fHistClusCent(0),
   fHistJetsCent(0),
@@ -110,18 +84,10 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
   fHistJetsParts(0),
   fHistCellsCent(0),
   fHistCellsTracks(0),
-  fHistMaxL1FastORCent(0),
-  fHistMaxL1ClusCent(0),
-  fHistMaxL1ThrCent(0),
-  fHistTracksPt(0),
-  fHistTrPhiEta(0),
   fHistTrEmcPhiEta(0),
   fHistTrPhiEtaNonProp(0),
   fHistDeltaEtaPt(0),
   fHistDeltaPhiPt(0),
-  fHistDeltaEtaNewProp(0),
-  fHistDeltaPhiNewProp(0),
-  fHistClusPhiEtaEnergy(0),
   fHistNCellsEnergy(0),
   fHistClusTimeEnergy(0),
   fHistCellsEnergy(0),
@@ -136,20 +102,14 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
     fHistTrackEta[i] = 0;
   }
 
-  for (Int_t i = 0; i < 6; i++) {
-    fHistTrackPhiPt[i] = 0;
-    fHistTrackEtaPt[i] = 0;
-  }
-
   for (Int_t i = 0; i < 4; i++) {
-    fHistJetsPhiEta[i] = 0;
-    fHistJetsPtNonBias[i] = 0;
-    fHistJetsPtTrack[i] = 0;
-    fHistJetsPtClus[i] = 0;
-    fHistJetsPt[i] = 0;
+    fHistTrPhiEtaPt[i] = 0;
+    fHistClusPhiEtaEnergy[i] = 0;
+    fHistJetsPhiEtaPt[i] = 0;
     fHistJetsPtArea[i] = 0;
-    fHistJetsPtAreaNonBias[i] = 0;
   }
+
+  SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
@@ -163,36 +123,73 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 {
   // Create histograms
 
-  OpenFile(1);
-  fOutput = new TList();
-  fOutput->SetOwner();
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 
-  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
-  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
-  fHistCentrality->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistCentrality);
+  if (!fTracksName.IsNull()) {
+    fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
+    fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
+    fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
+    fOutput->Add(fHistTracksCent);
 
-  fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
-  fHistZVertex->GetXaxis()->SetTitle("z");
-  fHistZVertex->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistZVertex);
+    TString histname;
 
-  fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
-  fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
-  fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
-  fOutput->Add(fHistTracksCent);
+    for (Int_t i = 0; i < 4; i++) {
+      histname = "fHistTrPhiEtaPt_";
+      histname += i;
+      fHistTrPhiEtaPt[i] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
+      fHistTrPhiEtaPt[i] ->GetXaxis()->SetTitle("#eta");
+      fHistTrPhiEtaPt[i] ->GetYaxis()->SetTitle("#phi");
+      fHistTrPhiEtaPt[i] ->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+      fOutput->Add(fHistTrPhiEtaPt[i]);
+    }
 
-  fHistJetsCent = new TH2F("fHistJetsCent","Jets vs. centrality", 100, 0, 100, 60, 0, 60);
-  fHistJetsCent->GetXaxis()->SetTitle("Centrality (%)");
-  fHistJetsCent->GetYaxis()->SetTitle("No. of jets");
-  fOutput->Add(fHistJetsCent);
+    fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+    fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
+    fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
+    fOutput->Add(fHistTrEmcPhiEta);
+
+    fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+    fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
+    fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
+    fOutput->Add(fHistTrPhiEtaNonProp);
+    
+    fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
+    fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
+    fOutput->Add(fHistDeltaEtaPt);
+    
+    fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
+    fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
+    fOutput->Add(fHistDeltaPhiPt);
+
+    for (Int_t i = 0; i < 5; i++) {
+      TString histnamephi("fHistTrackPhi_");
+      histnamephi += i;
+      fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 201, 0, TMath::Pi() * 2.01);
+      fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
+      fOutput->Add(fHistTrackPhi[i]);
+      
+      TString histnameeta("fHistTrackEta_");
+      histnameeta += i;
+      fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -1, 1);
+      fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
+      fOutput->Add(fHistTrackEta[i]);
+    }
 
-  fHistJetsParts = new TH2F("fHistJetsParts","Jets vs. centrality", fNbins, 0, 6000, 60, 0, 60);
-  fHistJetsParts->GetXaxis()->SetTitle("No. of particles");
-  fHistJetsParts->GetYaxis()->SetTitle("No. of jets");
-  fOutput->Add(fHistJetsParts);
+    fHistTrackPhi[0]->SetLineColor(kRed);
+    fHistTrackEta[0]->SetLineColor(kRed);
+    fHistTrackPhi[1]->SetLineColor(kBlue);
+    fHistTrackEta[1]->SetLineColor(kBlue);
+    fHistTrackPhi[2]->SetLineColor(kGreen);
+    fHistTrackEta[2]->SetLineColor(kGreen);
+    fHistTrackPhi[3]->SetLineColor(kOrange);
+    fHistTrackEta[3]->SetLineColor(kOrange);
+    fHistTrackPhi[4]->SetLineColor(kBlack);
+    fHistTrackEta[4]->SetLineColor(kBlack);
+  }
 
-  if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
+  if (!fCaloName.IsNull()) {
     fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
     fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
     fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
@@ -213,204 +210,78 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
     fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
     fOutput->Add(fHistCellsTracks);
 
+    TString histname;
+
+    for (Int_t i = 0; i < 4; i++) {
+      histname = "fHistClusPhiEtaEnergy_";
+      histname += i;
+      fHistClusPhiEtaEnergy[i] = new TH3F(histname, histname, 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
+      fHistClusPhiEtaEnergy[i]->GetXaxis()->SetTitle("#eta");
+      fHistClusPhiEtaEnergy[i]->GetYaxis()->SetTitle("#phi");
+      fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("Energy (GeV)");
+      fOutput->Add(fHistClusPhiEtaEnergy[i]);
+    }
+
     fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins,  -1e-6, 1e-6);
     fHistClusTimeEnergy->GetXaxis()->SetTitle("Energy (GeV)");
     fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
     fOutput->Add(fHistClusTimeEnergy);
 
-    if (fDoTrigger) {
-      fHistMaxL1FastORCent = new TH2F("fHistMaxL1FastORCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
-      fHistMaxL1FastORCent->GetXaxis()->SetTitle("Centrality [%]");
-      fHistMaxL1FastORCent->GetYaxis()->SetTitle("Maximum L1 FastOR");
-      fOutput->Add(fHistMaxL1FastORCent);
-      
-      fHistMaxL1ClusCent = new TH2F("fHistMaxL1ClusCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
-      fHistMaxL1ClusCent->GetXaxis()->SetTitle("Centrality [%]");
-      fHistMaxL1ClusCent->GetYaxis()->SetTitle("Maximum L1 trigger cluster");
-      fOutput->Add(fHistMaxL1ClusCent);
-      
-      fHistMaxL1ThrCent = new TH2F("fHistMaxL1ThrCent","fHistMaxL1ThrCent", 100, 0, 100, 250, 0, 250);
-      fHistMaxL1ThrCent->GetXaxis()->SetTitle("Centrality [%]");
-      fHistMaxL1ThrCent->GetYaxis()->SetTitle("Maximum L1 threshold");
-      fOutput->Add(fHistMaxL1ThrCent);
-    }
-  }
-
-  fHistTracksPt = new TH1F("fHistTracksPt","p_{T} spectrum of reconstructed tracks", fNbins, fMinBinPt, fMaxBinPt);
-  fHistTracksPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistTracksPt->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistTracksPt);
-       
-  fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
-  fHistTrPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistTrPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistTrPhiEta);
-
-  fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
-  fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistTrEmcPhiEta);
-
-  fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
-  fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
-  fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistTrPhiEtaNonProp);
-
-  fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
-  fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
-  fOutput->Add(fHistDeltaEtaPt);
-
-  fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
-  fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
-  fOutput->Add(fHistDeltaPhiPt);
-
-  if (fRepropagateTracks) {
-    fHistDeltaEtaNewProp = new TH1F("fHistDeltaEtaNewProp","fHistDeltaEtaNewProp", 800, -0.4, 0.4);
-    fHistDeltaEtaNewProp->GetXaxis()->SetTitle("#delta#eta");
-    fHistDeltaEtaNewProp->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaEtaNewProp);
-    
-    fHistDeltaPhiNewProp = new TH1F("fHistDeltaPhiNewProp","fHistDeltaPhiNewProp", 2560, -1.6, 3.2);
-    fHistDeltaPhiNewProp->GetXaxis()->SetTitle("#delta#phi");
-    fHistDeltaPhiNewProp->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPhiNewProp);
-  }
-
-  if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-    fHistClusPhiEtaEnergy = new TH3F("fHistClusPhiEtaEnergy","Phi-Eta-Energy distribution of clusters", 
-                                    fNbins, fMinBinPt, fMaxBinPt, 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01);
-    fHistClusPhiEtaEnergy->GetXaxis()->SetTitle("E [GeV]");
-    fHistClusPhiEtaEnergy->GetYaxis()->SetTitle("#eta");
-    fHistClusPhiEtaEnergy->GetZaxis()->SetTitle("#phi");
-    fOutput->Add(fHistClusPhiEtaEnergy);
-
     fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
-    fHistNCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
+    fHistNCellsEnergy->GetXaxis()->SetTitle("Energy (GeV)");
     fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
-    fOutput->Add(fHistNCellsEnergy);
-  }
-
-  if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-   
+    fOutput->Add(fHistNCellsEnergy); 
+     
     fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinBinPt, fMaxBinPt);
-    fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
+    fHistCellsEnergy->GetXaxis()->SetTitle("Energy (GeV)");
     fHistCellsEnergy->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistCellsEnergy);
     
     fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy", 
                                (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
-    fHistChVSneCells->GetXaxis()->SetTitle("E [GeV]");
-    fHistChVSneCells->GetYaxis()->SetTitle("P [GeV/c]");
+    fHistChVSneCells->GetXaxis()->SetTitle("Energy (GeV)");
+    fHistChVSneCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
     fOutput->Add(fHistChVSneCells);
     
     fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy", 
                               (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
-    fHistChVSneClus->GetXaxis()->SetTitle("E [GeV]");
-    fHistChVSneClus->GetYaxis()->SetTitle("P [GeV/c]");
+    fHistChVSneClus->GetXaxis()->SetTitle("Energy (GeV)");
+    fHistChVSneClus->GetYaxis()->SetTitle("Momentum (GeV/c)");
     fOutput->Add(fHistChVSneClus);
     
     fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy", 
                                    (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt , fMaxBinPt * 2.5);
-    fHistChVSneCorrCells->GetXaxis()->SetTitle("E [GeV]");
-    fHistChVSneCorrCells->GetYaxis()->SetTitle("P [GeV/c]");
+    fHistChVSneCorrCells->GetXaxis()->SetTitle("Energy (GeV)");
+    fHistChVSneCorrCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
     fOutput->Add(fHistChVSneCorrCells);
   }
-  for (Int_t i = 0; i < 5; i++) {
-    TString histnamephi("fHistTrackPhi_");
-    histnamephi += i;
-    fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 201, 0, TMath::Pi() * 2.01);
-    fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
-    fOutput->Add(fHistTrackPhi[i]);
-
-    TString histnameeta("fHistTrackEta_");
-    histnameeta += i;
-    fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -1, 1);
-    fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
-    fOutput->Add(fHistTrackEta[i]);
-  }
-
-  for (Int_t i = 0; i < 6; i++) {
-    TString histnamephipt("fHistTrackPhiPt_");
-    histnamephipt += i;
-    fHistTrackPhiPt[i] = new TH1F(histnamephipt.Data(),histnamephipt.Data(), 201, 0, TMath::Pi() * 2.01);
-    fHistTrackPhiPt[i]->GetXaxis()->SetTitle("Phi");
-    fOutput->Add(fHistTrackPhiPt[i]);
-
-    TString histnameetapt("fHistTrackEtaPt_");
-    histnameetapt += i;
-    fHistTrackEtaPt[i] = new TH1F(histnameetapt.Data(),histnameetapt.Data(), 100, -1, 1);
-    fHistTrackEtaPt[i]->GetXaxis()->SetTitle("Eta");
-    fOutput->Add(fHistTrackEtaPt[i]);
-  }
-  
-  fHistTrackPhi[0]->SetLineColor(kRed);
-  fHistTrackEta[0]->SetLineColor(kRed);
-  fHistTrackPhi[1]->SetLineColor(kBlue);
-  fHistTrackEta[1]->SetLineColor(kBlue);
-  fHistTrackPhi[2]->SetLineColor(kGreen);
-  fHistTrackEta[2]->SetLineColor(kGreen);
-  fHistTrackPhi[3]->SetLineColor(kOrange);
-  fHistTrackEta[3]->SetLineColor(kOrange);
-  fHistTrackPhi[4]->SetLineColor(kBlack);
-  fHistTrackEta[4]->SetLineColor(kBlack);
-
+       
   if (!fJetsName.IsNull()) {
+    fHistJetsCent = new TH2F("fHistJetsCent","Jets vs. centrality", 100, 0, 100, 150, -0.5, 149.5);
+    fHistJetsCent->GetXaxis()->SetTitle("Centrality (%)");
+    fHistJetsCent->GetYaxis()->SetTitle("No. of jets");
+    fOutput->Add(fHistJetsCent);
+    
+    fHistJetsParts = new TH2F("fHistJetsParts","Jets vs. centrality", fNbins, 0, 6000, 150, -0.5, 149.5);
+    fHistJetsParts->GetXaxis()->SetTitle("No. of particles");
+    fHistJetsParts->GetYaxis()->SetTitle("No. of jets");
+    fOutput->Add(fHistJetsParts);
 
     TString histname;
 
     for (Int_t i = 0; i < 4; i++) {
-      histname = "fHistJetsPhiEta_";
-      histname += i;
-      fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01);
-      fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
-      fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
-      fHistJetsPhiEta[i]->GetZaxis()->SetTitle("p_{T} [GeV/c]");
-      fOutput->Add(fHistJetsPhiEta[i]);
-
-      histname = "fHistJetsPtNonBias_";
-      histname += i;
-      fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
-      fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-      fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistJetsPtNonBias[i]);
-
-      histname = "fHistJetsPtTrack_";
-      histname += i;
-      fHistJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
-      fHistJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-      fHistJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistJetsPtTrack[i]);
-
-      if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-       histname = "fHistJetsPtClus_";
-       histname += i;
-       fHistJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
-       fHistJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-       fHistJetsPtClus[i]->GetYaxis()->SetTitle("counts");
-       fOutput->Add(fHistJetsPtClus[i]);
-      }
-
-      histname = "fHistJetsPt_";
-      histname += i;
-      fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
-      fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-      fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistJetsPt[i]);
-
-      histname = "fHistJetsPtAreaNonBias_";
+      histname = "fHistJetsPhiEtaPt_";
       histname += i;
-      fHistJetsPtAreaNonBias[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
-      fHistJetsPtAreaNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-      fHistJetsPtAreaNonBias[i]->GetYaxis()->SetTitle("area");
-      fOutput->Add(fHistJetsPtAreaNonBias[i]);
+      fHistJetsPhiEtaPt[i] = new TH3F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01,  (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
+      fHistJetsPhiEtaPt[i]->GetXaxis()->SetTitle("#eta");
+      fHistJetsPhiEtaPt[i]->GetYaxis()->SetTitle("#phi");
+      fHistJetsPhiEtaPt[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+      fOutput->Add(fHistJetsPhiEtaPt[i]);
 
       histname = "fHistJetsPtArea_";
       histname += i;
-      fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
-      fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+      fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
+      fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
       fOutput->Add(fHistJetsPtArea[i]);
     }
@@ -431,22 +302,6 @@ Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
   fNtracks = 0;
   fNjets = 0;
 
-  if (!fTrgClusName.IsNull() && fDoTrigger && !fTrgClusters) {
-    fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
-    if (!fTrgClusters) {
-      AliError(Form("%s: Could not retrieve trigger clusters %s!", GetName(), fTrgClusName.Data())); 
-      return kFALSE;
-    }
-    else {
-      TClass *cl = fTrgClusters->GetClass();
-      if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
-       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fTrgClusName.Data())); 
-       fTrgClusters = 0;
-       return kFALSE;
-      }
-    }
-  }
-
   return kTRUE;
 }
 
@@ -456,18 +311,17 @@ Bool_t AliAnalysisTaskSAQA::FillHistograms()
 {
   // Fill histograms.
 
-  fHistCentrality->Fill(fCent);
-  fHistZVertex->Fill(fVertex[2]);
+  Float_t clusSum = 0;
+  Float_t trackSum = 0;
 
-  Float_t trackSum = DoTrackLoop();
+  if (fTracks) {
+    trackSum = DoTrackLoop();
 
-  DoJetLoop();
-
-  fHistTracksCent->Fill(fCent, fNtracks);
-
-  if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
+    fHistTracksCent->Fill(fCent, fNtracks);
+  } 
 
-    Float_t clusSum = DoClusterLoop();
+  if (fCaloClusters) {
+    clusSum = DoClusterLoop();
 
     fHistClusCent->Fill(fCent, fNclusters);
     fHistClusTracks->Fill(fNtracks, fNclusters);
@@ -477,34 +331,22 @@ Bool_t AliAnalysisTaskSAQA::FillHistograms()
     Int_t ncells = DoCellLoop(cellSum, cellCutSum);
 
     if (fTracks)
-      fHistCellsTracks->Fill(fTracks->GetEntriesFast(), ncells);
+      fHistCellsTracks->Fill(fNtracks, ncells);
 
     fHistCellsCent->Fill(fCent, ncells);
     
     fHistChVSneCells->Fill(cellSum, trackSum);
     fHistChVSneClus->Fill(clusSum, trackSum);
     fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
+  }
 
-    if (fDoTrigger) {
-      Float_t maxTrgClus = DoTriggerClusLoop();
-      fHistMaxL1ClusCent->Fill(fCent, maxTrgClus);
+  if (fJets) {
+    DoJetLoop();
     
-      Int_t maxL1amp = -1;
-      Int_t maxL1thr = -1;
-    
-      DoTriggerPrimitives(maxL1amp, maxL1thr);
-      
-      if (maxL1amp > -1) 
-       fHistMaxL1FastORCent->Fill(fCent, maxL1amp);
-      
-      if (maxL1thr > -1) 
-       fHistMaxL1ThrCent->Fill(fCent, maxL1thr);
-    }
+    fHistJetsCent->Fill(fCent, fNjets);
+    fHistJetsParts->Fill(fNtracks + fNclusters, fNjets);
   }
 
-  fHistJetsCent->Fill(fCent, fNjets);
-  fHistJetsParts->Fill(fNtracks + fNclusters, fNjets);
-
   return kTRUE;
 }
 
@@ -560,7 +402,7 @@ Float_t AliAnalysisTaskSAQA::DoClusterLoop()
     TLorentzVector nPart;
     cluster->GetMomentum(nPart, fVertex);
 
-    fHistClusPhiEtaEnergy->Fill(cluster->E(), nPart.Eta(), nPart.Phi());
+    fHistClusPhiEtaEnergy[fCentBin]->Fill(nPart.Eta(), nPart.Phi(), cluster->E());
     fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
 
     fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
@@ -601,47 +443,20 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
       continue;
 
     fNtracks++;
-    
-    fHistTracksPt->Fill(track->Pt());
 
     sum += track->P();
+      
+    fHistTrPhiEtaPt[fCentBin]->Fill(track->Eta(), track->Phi(), track->Pt());
 
     Int_t label = track->GetLabel();
-      
-    fHistTrPhiEta->Fill(track->Eta(), track->Phi());
-    
-    fHistTrackEta[4]->Fill(track->Eta());
-    fHistTrackPhi[4]->Fill(track->Phi());
 
     if (label >= 0 && label < 4) {
       fHistTrackEta[label]->Fill(track->Eta());
       fHistTrackPhi[label]->Fill(track->Phi());
     }
 
-    if (track->Pt() < 0.5) {
-      fHistTrackPhiPt[0]->Fill(track->Phi());
-      fHistTrackEtaPt[0]->Fill(track->Eta());
-    }
-    else if (track->Pt() < 1) {
-      fHistTrackPhiPt[1]->Fill(track->Phi());
-      fHistTrackEtaPt[1]->Fill(track->Eta());
-    }
-    else if (track->Pt() < 2) {
-      fHistTrackPhiPt[2]->Fill(track->Phi());
-      fHistTrackEtaPt[2]->Fill(track->Eta());
-    }
-    else if (track->Pt() < 3) {
-      fHistTrackPhiPt[3]->Fill(track->Phi());
-      fHistTrackEtaPt[3]->Fill(track->Eta());
-    }
-    else if (track->Pt() < 5) {
-      fHistTrackPhiPt[4]->Fill(track->Phi());
-      fHistTrackEtaPt[4]->Fill(track->Eta());
-    }
-    else {
-      fHistTrackPhiPt[5]->Fill(track->Phi());
-      fHistTrackEtaPt[5]->Fill(track->Eta());
-    }
+    fHistTrackEta[4]->Fill(track->Eta());
+    fHistTrackPhi[4]->Fill(track->Phi());
 
     if (!vtrack)
       continue;
@@ -652,60 +467,11 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
     fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
     fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
     fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
-
-    if (fRepropagateTracks && vtrack->GetTrackEtaOnEMCal() > -2) {    
-      Float_t propeta = -999, propphi = -999;
-      PropagateTrack(vtrack, propeta, propphi);
-      fHistDeltaEtaNewProp->Fill(propeta - vtrack->GetTrackEtaOnEMCal());
-      fHistDeltaPhiNewProp->Fill(propphi - vtrack->GetTrackPhiOnEMCal());
-    }
   }
   
   return sum;
 }
 
-//____________________________________________________________________________
-void AliAnalysisTaskSAQA::PropagateTrack(AliVTrack *track, Float_t &eta, Float_t &phi)
-{
-  eta = -999;
-  phi = -999;
-
-  if (!track)
-    return;
-
-  if (track->Pt() == 0) 
-    return;
-
-  // init the magnetic field if not already on
-  if(!TGeoGlobalMagField::Instance()->GetField()) {
-    AliInfo("Init the magnetic field\n");
-    AliAODEvent* aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
-    if (aodevent) {
-      Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
-      Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
-      AliMagF *field  = AliMagF::CreateFieldMap(curSol,curDip);
-      TGeoGlobalMagField::Instance()->SetField(field);
-    }
-  }
-    
-  Double_t cv[21];
-  for (Int_t i = 0; i < 21; i++) cv[i] = 0;
-
-  Double_t pos[3], mom[3];
-  track->GetXYZ(pos);
-  track->GetPxPyPz(mom);
-  AliExternalTrackParam *trackParam = new AliExternalTrackParam(pos, mom, cv, track->Charge());
-  if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, 430., 0.139, 20, kTRUE, 0.8, -1)) return;
-  Double_t trkPos[3] = {0., 0., 0.};
-  if(!trackParam->GetXYZ(trkPos)) return;
-  TVector3 trkPosVec(trkPos[0], trkPos[1], trkPos[2]);
-  eta = trkPosVec.Eta();
-  phi = trkPosVec.Phi();
-  if(phi < 0)
-    phi += 2 * TMath::Pi();
-}
-
 //________________________________________________________________________
 void AliAnalysisTaskSAQA::DoJetLoop()
 {
@@ -725,83 +491,16 @@ void AliAnalysisTaskSAQA::DoJetLoop()
       continue;
     }  
 
-    if (!AcceptJet(jet, kFALSE))
+    if (!AcceptJet(jet))
       continue;
 
-    fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
-    fHistJetsPtAreaNonBias[fCentBin]->Fill(jet->Pt(), jet->Area());
-
     fNjets++;
 
-    if (jet->MaxTrackPt() > fPtBiasJetTrack)
-      fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
-    
-    if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus)
-      fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
-    
-    if (!AcceptBiasJet(jet))
-      continue;
-
-    fHistJetsPt[fCentBin]->Fill(jet->Pt());
+    fHistJetsPhiEtaPt[fCentBin]->Fill(jet->Eta(), jet->Phi(), jet->Pt());
     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
-
-    fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
   }
 }
 
-//________________________________________________________________________
-Float_t AliAnalysisTaskSAQA::DoTriggerClusLoop()
-{
-  // Do trigger cluster loop.
-
-  if (!fTrgClusters)
-    return 0;
-
-  Int_t ntrgclusters = fTrgClusters->GetEntriesFast();
-  Float_t maxTrgClus = 0;
-
-  for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
-    AliVCluster* cluster = static_cast<AliVCluster*>(fTrgClusters->At(iClusters));
-    if (!cluster) {
-      AliError(Form("Could not receive cluster %d", iClusters));
-      continue;
-    }  
-    
-    if (!(cluster->IsEMCAL())) continue;
-
-    if (cluster->E() > maxTrgClus)
-      maxTrgClus = cluster->E();
-
-  }
-  return maxTrgClus;
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSAQA::DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)
-{
-  // Do trigger primitives loop.
-
-  AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
-
-  if (!triggers || triggers->GetEntries() == 0)
-    return;
-    
-  triggers->Reset();
-  Int_t L1amp = 0;
-  Int_t L1thr = 0;
-  maxL1amp = -1;
-  maxL1thr = -1;
-
-  while (triggers->Next()) {
-    triggers->GetL1TimeSum(L1amp);
-    if (maxL1amp < L1amp) 
-      maxL1amp = L1amp;
-
-    triggers->GetL1Threshold(L1thr);
-    if (maxL1thr < L1thr) 
-      maxL1thr = L1thr;
-  }
-}
 
 //________________________________________________________________________
 void AliAnalysisTaskSAQA::Terminate(Option_t *) 
index 2afdfaafe3255f99febcac256627deab5280af32..a40766ff227a8d042af71488848d8d4c6c4ff5b7 100644 (file)
@@ -20,35 +20,23 @@ class AliAnalysisTaskSAQA : public AliAnalysisTaskEmcalJet {
   void                        UserCreateOutputObjects();
   void                        Terminate(Option_t *option);
 
-  void                        SetTrgClusName(const char *n)                        { fTrgClusName        = n          ; }
   void                        SetCellEnergyCut(Float_t cut)                        { fCellEnergyCut      = cut        ; }
-  void                        SetDoTrigger(Bool_t trg = kTRUE)                     { fDoTrigger          = trg        ; }
-  void                        SetDoRepropagateTracks(Bool_t p = kTRUE)             { fRepropagateTracks  = p          ; }
 
  protected:
 
   Bool_t                      FillHistograms()                                              ;
   Bool_t                      RetrieveEventObjects()                                        ;
   Int_t                       DoCellLoop(Float_t &sum, Float_t &sum_cut)                    ;
-  void                        DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)         ;
-  Float_t                     DoTriggerClusLoop()                                           ;
   Float_t                     DoTrackLoop()                                                 ;
   Float_t                     DoClusterLoop()                                               ;
   void                        DoJetLoop()                                                   ;
-  void                        PropagateTrack(AliVTrack *track, Float_t &eta, Float_t &phi)  ;
 
   Float_t                     fCellEnergyCut;            // Energy cell cut
-  Bool_t                      fDoTrigger;                // Make trigger qa plots
-  Bool_t                      fRepropagateTracks;        // Repropagate tracks to the EMCal surface
-  TString                     fTrgClusName;              // Name of trg clus name
-  TClonesArray               *fTrgClusters;              //!Trg Clusters
   Int_t                       fNclusters;                //!Number of accepted clusters in the event
   Int_t                       fNtracks;                  //!Number of accepted tracks in the event
   Int_t                       fNjets;                    //!Number of accepted jets in the event
  
   // General histograms
-  TH1F                       *fHistCentrality;           //!Event centrality distribution
-  TH1F                       *fHistZVertex;              //!Z vertex position
   TH2F                       *fHistTracksCent;           //!Number of tracks vs. centrality
   TH2F                       *fHistClusCent;             //!Number of clusters vs. centrality
   TH2F                       *fHistJetsCent;             //!Number of jets vs. centrality
@@ -56,47 +44,39 @@ class AliAnalysisTaskSAQA : public AliAnalysisTaskEmcalJet {
   TH2F                       *fHistJetsParts;            //!Number of jets vs. number of particles (tracks+clusters)
   TH2F                       *fHistCellsCent;            //!Number of cells vs. centrality
   TH2F                       *fHistCellsTracks;          //!Number of cells vs. number of tracks
-  // EMCAL trigger
-  TH2F                       *fHistMaxL1FastORCent;      //!Maximum L1 trigger FastOR amplitude vs. centrality
-  TH2F                       *fHistMaxL1ClusCent;        //!Maximum L1 trigger cluster amplitude vs. centrality
-  TH2F                       *fHistMaxL1ThrCent;         //!Maximum L1 trigger threshold vs. centrality
+
   // Tracks
-  TH1F                       *fHistTracksPt;             //!Pt spectrum of tracks
-  TH2F                       *fHistTrPhiEta;             //!Phi-Eta distribution of tracks
+  TH3F                       *fHistTrPhiEtaPt[4];        //!Phi-Eta-Pt distribution of tracks
   TH2F                       *fHistTrEmcPhiEta;          //!Phi-Eta emcal propagated distribution of tracks
   TH2F                       *fHistTrPhiEtaNonProp;      //!Phi-Eta distribution of non emcal propagated tracks
   TH2F                       *fHistDeltaEtaPt;           //!Eta-EtaProp vs. Pt
   TH2F                       *fHistDeltaPhiPt;           //!Phi-PhiProp vs. Pt
-  TH1F                       *fHistDeltaEtaNewProp;      //!NewEtaProp-EtaProp
-  TH1F                       *fHistDeltaPhiNewProp;      //!NewPhiProp-PhiProp
+
   // Clusters
-  TH3F                       *fHistClusPhiEtaEnergy;     //!Phi-Eta-Energy distribution of clusters
+  TH3F                       *fHistClusPhiEtaEnergy[4];  //!Phi-Eta-Energy distribution of clusters
   TH2F                       *fHistNCellsEnergy;         //!Number of cells vs. energy of cluster
   TH2F                       *fHistClusTimeEnergy;       //!Time vs. energy of cluster
+
   //Jets
-  TH2F                       *fHistJetsPhiEta[4];        //!Phi-Eta distribution of jets
-  TH1F                       *fHistJetsPtNonBias[4];     //!Non biased inclusive jet pt spectrum
-  TH1F                       *fHistJetsPtClus[4];        //!Inclusive jet pt spectrum cluster biased
-  TH1F                       *fHistJetsPtTrack[4];       //!Inclusive jet pt spectrum track biased
-  TH1F                       *fHistJetsPt[4];            //!Biased inclusive jet pt spectrum
-  TH2F                       *fHistJetsPtAreaNonBias[4]; //!Non biased pt vs. area of jets
-  TH2F                       *fHistJetsPtArea[4];        //!Biased pt vs. area of jets
+  TH3F                       *fHistJetsPhiEtaPt[4];      //!Phi-Eta distribution of jets
+  TH2F                       *fHistJetsPtArea[4];        //!Pt vs. area of jets
+
   // EMCAL Cells
   TH1F                       *fHistCellsEnergy;          //!Energy spectrum of cells
+
   // Had corr QA
   TH2F                       *fHistChVSneCells;          //!Charged vs. neutral (cells) energy
   TH2F                       *fHistChVSneClus;           //!Charged vs. neutral (clusters) energy
   TH2F                       *fHistChVSneCorrCells;      //!Charged vs. neutral (corrected cells) energy
+
   // Hybrid tracks
   TH1F                       *fHistTrackPhi[5];          //!Phi distribution of hybrid tracks
   TH1F                       *fHistTrackEta[5];          //!Eta distribution of hybrid tracks
-  TH1F                       *fHistTrackPhiPt[6];        //!Phi distribution of hybrid tracks per pt bins
-  TH1F                       *fHistTrackEtaPt[6];        //!Eta distribution of hybrid tracks per pt bins
 
  private:
   AliAnalysisTaskSAQA(const AliAnalysisTaskSAQA&);            // not implemented
   AliAnalysisTaskSAQA &operator=(const AliAnalysisTaskSAQA&); // not implemented
 
-  ClassDef(AliAnalysisTaskSAQA, 13) // Quality task for Emcal analysis
+  ClassDef(AliAnalysisTaskSAQA, 14) // Quality task for Emcal analysis
 };
 #endif
diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskDeltaPt.C b/PWGJE/EMCALJetTasks/macros/AddTaskDeltaPt.C
new file mode 100644 (file)
index 0000000..3795788
--- /dev/null
@@ -0,0 +1,81 @@
+
+AliAnalysisTaskDeltaPt* AddTaskDeltaPt(
+  const char *ntracks            = "Tracks",
+  const char *nclusters          = "CaloClusters",
+  const char *njets              = "Jets",
+  const char *nembtracks         = "TracksEmbedded",
+  const char *nembclusters       = "CaloClustersEmbedded",
+  const char *nembjets           = "EmbJets",
+  const char *nrandtracks        = "TracksRandomized",
+  const char *nrandclusters      = "CaloClustersRandomized",
+  const char *nrho               = "Rho",
+  Double_t    jetradius          = 0.2,
+  Double_t    jetptcut           = 1,
+  Double_t    jetareacut         = 0.8,
+  Double_t    ptcut              = 0.15,
+  UInt_t      type               = AliAnalysisTaskEmcal::kTPC,
+  const char *taskname           = "AliAnalysisTaskDeltaPt"
+)
+{  
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+  {
+    ::Error("AddTaskSAJF", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler())
+  {
+    ::Error("AddTaskSAJF", "This task requires an input event handler");
+    return NULL;
+  }
+  
+  //-------------------------------------------------------
+  // Init the task and do settings
+  //-------------------------------------------------------
+  TString name(Form("%s_%s_%s_%s_R0%d_",taskname,ntracks,nclusters,nrho,(Int_t)floor(jetradius*100+0.5)));
+  if (type == AliAnalysisTaskEmcal::kTPC) 
+    name += "TPC";
+  else if (type == AliAnalysisTaskEmcal::kEMCAL) 
+    name += "EMCAL";
+  else if (type == AliAnalysisTaskEmcal::kUser) 
+    name += "USER";
+
+  AliAnalysisTaskDeltaPt* jetTask = new AliAnalysisTaskDeltaPt(name);
+  jetTask->SetAnaType(type);
+  jetTask->SetTracksName(ntracks);
+  jetTask->SetClusName(nclusters);
+  jetTask->SetJetsName(njets);
+  jetTask->SetEmbTracksName(nembtracks);
+  jetTask->SetEmbClusName(nembclusters);
+  jetTask->SetEmbJetsName(nembjets);
+  jetTask->SetRandTracksName(nrandtracks);
+  jetTask->SetRandClusName(nrandclusters);
+  jetTask->SetRhoName(nrho);
+  jetTask->SetPtCut(ptcut);
+  jetTask->SetJetRadius(jetradius);
+  jetTask->SetJetPtCut(jetptcut);
+  jetTask->SetPercAreaCut(jetareacut);
+  
+  //-------------------------------------------------------
+  // Final settings, pass to manager and set the containers
+  //-------------------------------------------------------
+  
+  mgr->AddTask(jetTask);
+  
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
+  TString contname(name);
+  contname += "_histos";
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(), 
+                                                           TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                           Form("%s", AliAnalysisManager::GetCommonFileName()));
+  mgr->ConnectInput  (jetTask, 0,  cinput1 );
+  mgr->ConnectOutput (jetTask, 1, coutput1 );
+  
+  return jetTask;
+}
index ae8e01b8bb1b6f6313442624e90ffb43c9459a41..d6d44466336630e1d2e2c881205c84b2dbf437a1 100644 (file)
@@ -1,18 +1,19 @@
 // $Id$
 
 AliAnalysisTaskRho* AddTaskRho(
-   const char *outfilename    = "AnalysisResults.root",
-   const char *nJets          = "Jets",
-   const char *nTracks        = "PicoTracks",   
-   const char *nRho           = "Rho",
-   const Double_t minPhi      = 0,
-   const Double_t maxPhi      = 2 * TMath::Pi(),
-   const Double_t minEta      = -0.3,
-   const Double_t maxEta      = 0.3,
-   const Double_t minArea     = 0.01,
-   const UInt_t   exclJets    = 1,
+   const char    *nJets       = "Jets",
+   const char    *nTracks     = "PicoTracks",
+   const char    *nClusters   = "CaloClusters",  
+   const char    *nRho        = "Rho",
+   Double_t       jetradius   = 0.2,
+   UInt_t         type        = AliAnalysisTaskEmcal::kTPC,
+   Double_t       jetareacut  = 0.01,
+   Double_t       emcareacut  = 0,
+   Double_t       ptcut       = 0.15,
+   TF1           *sfunc       = 0,
+   const UInt_t   exclJets    = 2,
    const Bool_t   histo       = kFALSE,
-   const char *taskname       = "Rho"
+   const char    *taskname    = "Rho"
 )
 {  
   // Get the pointer to the existing analysis manager via the static access method.
@@ -36,14 +37,25 @@ AliAnalysisTaskRho* AddTaskRho(
   // Init the task and do settings
   //-------------------------------------------------------
 
-  TString name(Form("%s_%s", taskname, nJets));
+  TString name(Form("%s_%s_", taskname, nJets));
+  if (type == AliAnalysisTaskEmcal::kTPC) 
+    name += "TPC";
+  else if (type == AliAnalysisTaskEmcal::kEMCAL) 
+    name += "EMCAL";
+  else if (type == AliAnalysisTaskEmcal::kUser) 
+    name += "USER";
   AliAnalysisTaskRho *rhotask = new AliAnalysisTaskRho(name, histo);
+  rhotask->SetAnaType(type);
+  rhotask->SetScaleFunction(sfunc);
   rhotask->SetJetsName(nJets);
   rhotask->SetTracksName(nTracks);
+  rhotask->SetClusName(nClusters);
   rhotask->SetRhoName(nRho);
-  rhotask->SetJetPhi(minPhi,maxPhi);
-  rhotask->SetJetEta(minEta,maxEta);
-  rhotask->SetAreaCut(minArea);
+  rhotask->SetJetAreaCut(jetareacut);
+  rhotask->SetAreaEmcCut(emcareacut);
+  rhotask->SetPtCut(ptcut);
+  rhotask->SetJetPtCut(0);
+  rhotask->SetJetRadius(jetradius);
   rhotask->SetExcludeLeadJets(exclJets);
 
   //-------------------------------------------------------
@@ -53,13 +65,14 @@ AliAnalysisTaskRho* AddTaskRho(
   mgr->AddTask(rhotask);
 
   // Create containers for input/output
-  mgr->ConnectInput (rhotask, 0, mgr->GetCommonInputContainer() );
+  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
   if (histo) {
-    AliAnalysisDataContainer *corho = mgr->CreateContainer(name,
-                                                           TList::Class(),
-                                                           AliAnalysisManager::kOutputContainer,
-                                                           outfilename);
-    mgr->ConnectOutput(rhotask, 1, corho);
+    TString contname(name);
+    contname += "_histos";
+    AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(), 
+                                                             TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                             Form("%s", AliAnalysisManager::GetCommonFileName()));
+    mgr->ConnectOutput(rhotask, 1, coutput1);
   }
 
   return rhotask;
index af8601d1c0551bc9bd9d0ce5ebaaae9922c2cbf8..ae9f3b4b433d86eeb17508cebec1665b6980eee0 100644 (file)
@@ -1,16 +1,20 @@
 // $Id$
 
 AliAnalysisTaskRhoAverage* AddTaskRhoAverage(
-   const char *nJets          = "Jets",
-   const char *nTracks        = "PicoTracks",   
-   const char *nClusters      = "CaloClustersCorr",  
-   const char *nRho           = "RhoAverage",
-   const Double_t minPhi      = 0,
-   const Double_t maxPhi      = 2 * TMath::Pi(),
-   const Double_t minEta      = -0.9,
-   const Double_t maxEta      = 0.9,
-   const Double_t minPt       = 0.15,
-   const char *taskname       = "RhoAverage"
+   const char    *nJets       = "Jets",
+   const char    *nTracks     = "PicoTracks",
+   const char    *nClusters   = "CaloClusters",  
+   const char    *nRho        = "Rho",
+   Double_t       jetradius   = 0.2,
+   UInt_t         type        = AliAnalysisTaskEmcal::kTPC,
+   Double_t       jetareacut  = 0.01,
+   Double_t       emcareacut  = 0,
+   Double_t       ptcut       = 0.15,
+   TF1           *sfunc       = 0,
+   const UInt_t   exclPart    = 2,
+   const UInt_t   rhotype     = 1,
+   const Bool_t   histo       = kFALSE,
+   const char    *taskname    = "RhoAverage"
 )
 {  
   // Get the pointer to the existing analysis manager via the static access method.
@@ -34,15 +38,27 @@ AliAnalysisTaskRhoAverage* AddTaskRhoAverage(
   // Init the task and do settings
   //-------------------------------------------------------
 
-  TString name(Form("%s_%s_%s", taskname, nTracks, nClusters));
-  AliAnalysisTaskRhoAverage *rhotask = new AliAnalysisTaskRhoAverage(name);
+  TString name(Form("%s_%s_%s_", taskname, nTracks, nClusters));
+  if (type == AliAnalysisTaskEmcal::kTPC) 
+    name += "TPC";
+  else if (type == AliAnalysisTaskEmcal::kEMCAL) 
+    name += "EMCAL";
+  else if (type == AliAnalysisTaskEmcal::kUser) 
+    name += "USER";
+  AliAnalysisTaskRhoAverage *rhotask = new AliAnalysisTaskRhoAverage(name, histo);
+  rhotask->SetAnaType(type);
+  rhotask->SetScaleFunction(sfunc);
   rhotask->SetJetsName(nJets);
   rhotask->SetTracksName(nTracks);
-  rhotask->SetClustersName(nClusters);
+  rhotask->SetClusName(nClusters);
   rhotask->SetRhoName(nRho);
-  rhotask->SetPhiLimits(minPhi,maxPhi);
-  rhotask->SetEtaLimits(minEta,maxEta);
-  rhotask->SetPtMin(minPt);
+  rhotask->SetJetAreaCut(jetareacut);
+  rhotask->SetAreaEmcCut(emcareacut);
+  rhotask->SetPtCut(ptcut);
+  rhotask->SetJetPtCut(0);
+  rhotask->SetJetRadius(jetradius);
+  rhotask->SetExcludeLeadPart(exclPart);
+  rhotask->SetRhoType(rhotype);
 
   //-------------------------------------------------------
   // Final settings, pass to manager and set the containers
@@ -51,7 +67,16 @@ AliAnalysisTaskRhoAverage* AddTaskRhoAverage(
   mgr->AddTask(rhotask);
 
   // Create containers for input/output
-  mgr->ConnectInput (rhotask, 0, mgr->GetCommonInputContainer() );
+  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
+
+  if (histo) {
+    TString contname(name);
+    contname += "_histos";
+    AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(), 
+                                                             TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                             Form("%s", AliAnalysisManager::GetCommonFileName()));
+    mgr->ConnectOutput(rhotask, 1, coutput1);
+  }
 
   return rhotask;
 }
index f5c3203f4ef6dcc617f4fef3ebe3918484065b2f..0b7d26b96aa17e8b6797c74054914c0e25787398 100644 (file)
@@ -1,9 +1,20 @@
 // $Id$
 
 AliAnalysisTaskRhoBase* AddTaskRhoBase(
-   const char *rhoname        = "Rho",
-   TF1        *rfunc          = 0,
-   const char *taskname       = "Rho_Base"
+   const char    *nJets       = "Jets",
+   const char    *nTracks     = "PicoTracks",
+   const char    *nClusters   = "CaloClusters",  
+   const char    *nRho        = "Rho",
+   Double_t       jetradius   = 0.2,
+   UInt_t         type        = AliAnalysisTaskEmcal::kTPC,
+   Double_t       jetareacut  = 0.01,
+   Double_t       emcareacut  = 0,
+   Double_t       ptcut       = 0.15,
+   TF1           *sfunc       = 0,
+   TF1           *rfunc       = 0,
+   const Bool_t   histo       = kFALSE,
+   const char    *taskname    = "RhoBase"
+
 )
 {  
   // Get the pointer to the existing analysis manager via the static access method.
@@ -27,9 +38,19 @@ AliAnalysisTaskRhoBase* AddTaskRhoBase(
   // Init the task and do settings
   //-------------------------------------------------------
 
-  AliAnalysisTaskRhoBase *rhotask = new AliAnalysisTaskRhoBase(taskname);
-  rhotask->SetRhoName(rhoname);
+  AliAnalysisTaskRhoBase *rhotask = new AliAnalysisTaskRhoBase(taskname,histo);
+  rhotask->SetAnaType(type);
   rhotask->SetRhoFunction(rfunc);
+  rhotask->SetScaleFunction(sfunc);
+  rhotask->SetJetsName(nJets);
+  rhotask->SetTracksName(nTracks);
+  rhotask->SetClusName(nClusters);
+  rhotask->SetRhoName(nRho);
+  rhotask->SetJetAreaCut(jetareacut);
+  rhotask->SetAreaEmcCut(emcareacut);
+  rhotask->SetPtCut(ptcut);
+  rhotask->SetJetPtCut(0);
+  rhotask->SetJetRadius(jetradius);
 
   //-------------------------------------------------------
   // Final settings, pass to manager and set the containers
@@ -38,7 +59,16 @@ AliAnalysisTaskRhoBase* AddTaskRhoBase(
   mgr->AddTask(rhotask);
 
   // Create containers for input/output
-  mgr->ConnectInput (rhotask, 0, mgr->GetCommonInputContainer() );
+  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
+
+  if (histo) {
+    TString contname(taskname);
+    contname += "_histos";
+    AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(), 
+                                                             TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                             Form("%s", AliAnalysisManager::GetCommonFileName()));
+    mgr->ConnectOutput(rhotask, 1, coutput1);
+  }
 
   return rhotask;
 }
index bf86424f3bb2a5359218000a36e86710289a65c2..631a565a1826e90988b0f3548f08573f8dfacb4a 100644 (file)
@@ -4,19 +4,13 @@ AliAnalysisTaskSAJF* AddTaskSAJF(
   const char *ntracks            = "Tracks",
   const char *nclusters          = "CaloClusters",
   const char *njets              = "Jets",
-  const char *nembtracks         = "TracksEmbedded",
-  const char *nembclusters       = "CaloClustersEmbedded",
-  const char *nembjets           = "EmbJets",
-  const char *nrandtracks        = "TracksRandomized",
-  const char *nrandclusters      = "CaloClustersRandomized",
   const char *nrho               = "Rho",
-  Double_t    jetradius          = 0.4,
+  Double_t    jetradius          = 0.2,
   Double_t    jetptcut           = 1,
   Double_t    jetareacut         = 0.8,
   Double_t    ptcut              = 0.15,
-  Double_t    jetBiasTrack       = 5,
-  Double_t    jetBiasClus        = 5,
   UInt_t      type               = AliAnalysisTaskEmcal::kTPC,
+  Int_t       leadhadtype        = 0,
   const char *taskname           = "AliAnalysisTaskSAJF"
 )
 {  
@@ -41,45 +35,24 @@ AliAnalysisTaskSAJF* AddTaskSAJF(
   // Init the task and do settings
   //-------------------------------------------------------
 
-  TString name(taskname);
-  name += "_";
-  name += njets;
-  name += "_";
-  name += nrho;
-  name += "_Track";
-  name += jetBiasTrack;
-  name += "_Clus";
-  name += jetBiasClus;
-  name += "_R0";
-  name += floor(jetradius*10+0.5);
-  name += "_PtCut";
-  name += floor(ptcut*1000+0.5);
-  name += "_";
+  TString name(Form("%s_%s_%s_R0%d_",taskname,njets,nrho,(Int_t)floor(jetradius*100+0.5)));
   if (type == AliAnalysisTaskEmcal::kTPC) 
     name += "TPC";
   else if (type == AliAnalysisTaskEmcal::kEMCAL) 
     name += "EMCAL";
-  else if (type == AliAnalysisTaskEmcal::kTPCSmall) 
-    name += "TPCSmall";
-  else if (type == AliAnalysisTaskEmcal::kEMCALOnly) 
-    name += "EMCALOnly";
+  else if (type == AliAnalysisTaskEmcal::kUser) 
+    name += "USER";
   AliAnalysisTaskSAJF* jetTask = new AliAnalysisTaskSAJF(name);
   jetTask->SetAnaType(type);
   jetTask->SetTracksName(ntracks);
   jetTask->SetClusName(nclusters);
   jetTask->SetJetsName(njets);
-  jetTask->SetEmbTracksName(nembtracks);
-  jetTask->SetEmbClusName(nembclusters);
-  jetTask->SetEmbJetsName(nembjets);
-  jetTask->SetRandTracksName(nrandtracks);
-  jetTask->SetRandClusName(nrandclusters);
   jetTask->SetRhoName(nrho);
   jetTask->SetPtCut(ptcut);
   jetTask->SetJetRadius(jetradius);
   jetTask->SetJetPtCut(jetptcut);
   jetTask->SetPercAreaCut(jetareacut);
-  jetTask->SetPtBiasJetTrack(jetBiasTrack);
-  jetTask->SetPtBiasJetClus(jetBiasClus);
+  jetTask->SetLeadingHadronType(leadhadtype);
   
   //-------------------------------------------------------
   // Final settings, pass to manager and set the containers
index 00b6d609f2a1d7fcdb8b109fbd2bed35db0add72..98eb07037b5f316d14b0dc038266a9952aa2936c 100644 (file)
@@ -4,13 +4,10 @@ AliAnalysisTaskSAQA* AddTaskSAQA(
   const char *ntracks            = "Tracks",
   const char *nclusters          = "CaloClusters",
   const char *njets              = "Jets",
-  const char *ntrgclusters       = "",
-  Double_t    jetradius          = 0.4,
+  Double_t    jetradius          = 0.2,
   Double_t    jetptcut           = 1,
   Double_t    jetareacut         = 0.8,
   Double_t    ptcut              = 0.15,
-  Double_t    jetBiasTrack       = 5,
-  Double_t    jetBiasClus        = 5,
   UInt_t      type               = AliAnalysisTaskEmcal::kTPC,
   const char *taskname           = "AliAnalysisTaskSAQA"
 )
@@ -35,37 +32,21 @@ AliAnalysisTaskSAQA* AddTaskSAQA(
   //-------------------------------------------------------
   // Init the task and do settings
   //-------------------------------------------------------
-  TString name(taskname);
-  name += "_";
-  name += njets;
-  name += "_Track";
-  name += jetBiasTrack;
-  name += "_Clus";
-  name += jetBiasClus;
-  name += "_R0";
-  name += floor(jetradius*10+0.5);
-  name += "_PtCut";
-  name += floor(ptcut*1000+0.5);
-  name += "_";
+  TString name(Form("%s_%s_%s_%s_R0%d_",taskname,njets,ntracks,nclusters,(Int_t)floor(jetradius*100+0.5)));
   if (type == AliAnalysisTaskEmcal::kTPC) 
     name += "TPC";
   else if (type == AliAnalysisTaskEmcal::kEMCAL) 
     name += "EMCAL";
-  else if (type == AliAnalysisTaskEmcal::kTPCSmall) 
-    name += "TPCSmall";
-  else if (type == AliAnalysisTaskEmcal::kEMCALOnly) 
-    name += "EMCALOnly";
+  else if (type == AliAnalysisTaskEmcal::kUser) 
+    name += "USER";
   AliAnalysisTaskSAQA* qaTask = new AliAnalysisTaskSAQA(name);
   qaTask->SetTracksName(ntracks);
   qaTask->SetClusName(nclusters);
   qaTask->SetJetsName(njets);
-  qaTask->SetTrgClusName(ntrgclusters);
   qaTask->SetJetRadius(jetradius);
   qaTask->SetJetPtCut(jetptcut);
   qaTask->SetPercAreaCut(jetareacut);
   qaTask->SetPtCut(ptcut);
-  qaTask->SetPtBiasJetTrack(jetBiasTrack);
-  qaTask->SetPtBiasJetClus(jetBiasClus);
   qaTask->SetAnaType(type);
 
   //-------------------------------------------------------
index 1152905e8036531eeb1b79b13c607303a7eb45ba..f101f0a104bf81b06b1fb83939356ea82f944d76 100644 (file)
@@ -7,9 +7,10 @@
 #pragma link off all functions;
 
 #pragma link C++ class AliAnalysisTaskEmcalJet+;
+#pragma link C++ class AliAnalysisTaskRhoBase+;
 #pragma link C++ class AliAnalysisTaskRho+;
 #pragma link C++ class AliAnalysisTaskRhoAverage+;
-#pragma link C++ class AliAnalysisTaskRhoBase+;
+#pragma link C++ class AliAnalysisTaskDeltaPt+;
 #pragma link C++ class AliAnalysisTaskScale+;
 #pragma link C++ class AliEmcalJet+;
 #pragma link C++ class AliEmcalJetTask+;