]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx
temporary patch to catch undermined runnumbers
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskSAQA.cxx
index c3160eba4538ea9988d0b16ff78d6756580e77ce..3353ab6bd6e28554e961831d16d407eb27413de1 100644 (file)
@@ -4,11 +4,11 @@
 //
 // Author: S.Aiola
 
-#include <TChain.h>
 #include <TClonesArray.h>
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TH3F.h>
+#include <THnSparse.h>
 #include <TList.h>
 #include <TLorentzVector.h>
 
 #include "AliExternalTrackParam.h"
 #include "AliTrackerBase.h"
 #include "AliLog.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALGeoParams.h"
+#include "AliPicoTrack.h"
+#include "AliVVZERO.h"
+#include "AliESDUtils.h"
 
 #include "AliAnalysisTaskSAQA.h"
 
@@ -32,104 +37,104 @@ ClassImp(AliAnalysisTaskSAQA)
 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() : 
   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
   fCellEnergyCut(0.1),
-  fDoTrigger(kFALSE),
-  fRepropagateTracks(kFALSE),
-  fTrgClusName("ClustersL1GAMMAFEE"),
-  fTrgClusters(0),
-  fHistCentrality(0),
-  fHistZVertex(0),
-  fHistTracksCent(0),
-  fHistClusCent(0),
-  fHistClusTracks(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),
+  fParticleLevel(kFALSE),
+  fIsMC(kFALSE),
+  fCentMethod2(""),
+  fCentMethod3(""),
+  fDoV0QA(0),
+  fDoEPQA(0),
+  fDoLeadingObjectPosition(0),
+  fMaxCellsInCluster(30),
+  fCent2(0),
+  fCent3(0),
+  fVZERO(0),
+  fV0ATotMult(0),
+  fV0CTotMult(0),
+  fHistEventQA(0),
+  fHistTrNegativeLabels(0),
+  fHistTrZeroLabels(0),
   fHistNCellsEnergy(0),
+  fHistFcrossEnergy(0),
   fHistClusTimeEnergy(0),
-  fHistCellsEnergy(0),
+  fHistClusEnergyMinusCellEnergy(0),
+  fHistCellsAbsIdEnergy(0),
   fHistChVSneCells(0),
   fHistChVSneClus(0),
   fHistChVSneCorrCells(0)
 {
   // Default constructor.
 
-  for (Int_t i = 0; i < 5; i++) {
-    fHistTrackPhi[i] = 0;
-    fHistTrackEta[i] = 0;
-  }
-
   for (Int_t i = 0; i < 4; i++) {
+    for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+    fHistTrPhiEtaZeroLab[i] = 0;
+    fHistTrPtZeroLab[i] = 0;
+    fHistTrEmcPhiEta[i] = 0;
+    fHistTrEmcPt[i] = 0;
+    fHistTrPhiEtaNonProp[i] = 0;
+    fHistTrPtNonProp[i] = 0;
+    fHistDeltaEtaPt[i] = 0;
+    fHistDeltaPhiPt[i] = 0;
+    fHistDeltaPtvsPt[i] = 0;
+    fHistClusPhiEtaEnergy[i] = 0;
+    fHistClusDeltaPhiEPEnergy[i] = 0;
+    fHistClusMCEnergyFraction[i] = 0;
     fHistJetsPhiEta[i] = 0;
-    fHistJetsPtNonBias[i] = 0;
-    fHistJetsPtTrack[i] = 0;
-    fHistJetsPtClus[i] = 0;
-    fHistJetsPt[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),
-  fHistCentrality(0),
-  fHistZVertex(0),
-  fHistTracksCent(0),
-  fHistClusCent(0),
-  fHistClusTracks(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),
+  fParticleLevel(kFALSE),
+  fIsMC(kFALSE),
+  fCentMethod2(""),
+  fCentMethod3(""),
+  fDoV0QA(0),
+  fDoEPQA(0),
+  fDoLeadingObjectPosition(0),
+  fMaxCellsInCluster(30),
+  fCent2(0),
+  fCent3(0),
+  fVZERO(0),
+  fV0ATotMult(0),
+  fV0CTotMult(0),
+  fHistEventQA(0),
+  fHistTrNegativeLabels(0),
+  fHistTrZeroLabels(0),
   fHistNCellsEnergy(0),
+  fHistFcrossEnergy(0),
   fHistClusTimeEnergy(0),
-  fHistCellsEnergy(0),
+  fHistClusEnergyMinusCellEnergy(0),
+  fHistCellsAbsIdEnergy(0),
   fHistChVSneCells(0),
   fHistChVSneClus(0),
   fHistChVSneCorrCells(0)
 {
   // Standard constructor.
 
-  for (Int_t i = 0; i < 5; i++) {
-    fHistTrackPhi[i] = 0;
-    fHistTrackEta[i] = 0;
-  }
-
   for (Int_t i = 0; i < 4; i++) {
+    for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+    fHistTrPhiEtaZeroLab[i] = 0;
+    fHistTrPtZeroLab[i] = 0;
+    fHistTrEmcPhiEta[i] = 0;
+    fHistTrEmcPt[i] = 0;
+    fHistTrPhiEtaNonProp[i] = 0;
+    fHistTrPtNonProp[i] = 0;
+    fHistDeltaEtaPt[i] = 0;
+    fHistDeltaPhiPt[i] = 0;
+    fHistDeltaPtvsPt[i] = 0;
+    fHistClusPhiEtaEnergy[i] = 0;
+    fHistClusDeltaPhiEPEnergy[i] = 0;
+    fHistClusMCEnergyFraction[i] = 0;
     fHistJetsPhiEta[i] = 0;
-    fHistJetsPtNonBias[i] = 0;
-    fHistJetsPtTrack[i] = 0;
-    fHistJetsPtClus[i] = 0;
-    fHistJetsPt[i] = 0;
     fHistJetsPtArea[i] = 0;
-    fHistJetsPtAreaNonBias[i] = 0;
   }
+
+  SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
@@ -143,237 +148,378 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 {
   // Create histograms
 
-  OpenFile(1);
-  fOutput = new TList();
-  fOutput->SetOwner();
-
-  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
-  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
-  fHistCentrality->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistCentrality);
-
-  fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
-  fHistZVertex->GetXaxis()->SetTitle("z");
-  fHistZVertex->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistZVertex);
-
-  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);
-
-  if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-    fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
-    fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
-    fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
-    fOutput->Add(fHistClusCent);
-
-    fHistClusTracks = new TH2F("fHistClusTracks","Clusters vs. tracks", fNbins, 0, 4000, fNbins, 0, 2000);
-    fHistClusTracks->GetXaxis()->SetTitle("No. of tracks");
-    fHistClusTracks->GetYaxis()->SetTitle("No. of clusters");
-    fOutput->Add(fHistClusTracks);
-
-    fHistCellsCent = new TH2F("fHistCellsCent","Cells vs. centrality", 100, 0, 100, fNbins, 0, 6000);
-    fHistCellsCent->GetXaxis()->SetTitle("Centrality (%)");
-    fHistCellsCent->GetYaxis()->SetTitle("No. of EMCal cells");
-    fOutput->Add(fHistCellsCent);
-
-    fHistCellsTracks = new TH2F("fHistCellsTracks","Cells vs. tracks", fNbins, 0, 4000, fNbins, 0, 6000);
-    fHistCellsTracks->GetXaxis()->SetTitle("No. of tracks");
-    fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
-    fOutput->Add(fHistCellsTracks);
-
-    fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins,  0, 1e-6);
-    fHistClusTimeEnergy->GetXaxis()->SetTitle("Energy (GeV)");
-    fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
-    fOutput->Add(fHistClusTimeEnergy);
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+  if (fParticleCollArray.GetEntriesFast()>0) {
+    if (!fParticleLevel && fIsMC) {
+      fHistTrNegativeLabels = new TH1F("fHistTrNegativeLabels","fHistTrNegativeLabels", 500, 0, 1);
+      fHistTrNegativeLabels->GetXaxis()->SetTitle("% of negative labels");
+      fHistTrNegativeLabels->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistTrNegativeLabels);
 
-    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);
+      fHistTrZeroLabels = new TH1F("fHistTrZeroLabels","fHistTrZeroLabels", 500, 0, 1);
+      fHistTrZeroLabels->GetXaxis()->SetTitle("% of negative labels");
+      fHistTrZeroLabels->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistTrZeroLabels);
     }
-  }
 
-  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", 80, -2, 2, 128, 0, 6.4);
-  fHistTrPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistTrPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistTrPhiEta);
-
-  fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 80, -2, 2, 128, 0, 6.4);
-  fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistTrEmcPhiEta);
-
-  fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 80, -2, 2, 128, 0, 6.4);
-  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);
+    TString histname;
+
+    Int_t nlabels = 4;
+    if (fParticleLevel)
+      nlabels = 1;
+
+    for (Int_t i = 0; i < fNcentBins; i++) {
+      for (Int_t j = 0; j < nlabels; j++) {
+       histname = Form("fHistTrPhiEtaPt_%d_%d",i,j);
+       fHistTrPhiEtaPt[i][j] = new TH3F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02, fNbins, fMinBinPt, fMaxBinPt);
+       fHistTrPhiEtaPt[i][j]->GetXaxis()->SetTitle("#eta");
+       fHistTrPhiEtaPt[i][j]->GetYaxis()->SetTitle("#phi");
+       fHistTrPhiEtaPt[i][j]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+       fOutput->Add(fHistTrPhiEtaPt[i][j]);
+      }
+
+      if (!fParticleLevel) {
+       if (fIsMC) {
+         histname = Form("fHistTrPhiEtaZeroLab_%d",i);
+         fHistTrPhiEtaZeroLab[i] = new TH2F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02);
+         fHistTrPhiEtaZeroLab[i]->GetXaxis()->SetTitle("#eta");
+         fHistTrPhiEtaZeroLab[i]->GetYaxis()->SetTitle("#phi");
+         fHistTrPhiEtaZeroLab[i]->GetZaxis()->SetTitle("counts");
+         fOutput->Add(fHistTrPhiEtaZeroLab[i]);
+
+         histname = Form("fHistTrPtZeroLab_%d",i);
+         fHistTrPtZeroLab[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
+         fHistTrPtZeroLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+         fHistTrPtZeroLab[i]->GetYaxis()->SetTitle("counts");
+         fOutput->Add(fHistTrPtZeroLab[i]);
+       }
+       
+       histname = Form("fHistTrEmcPhiEta_%d",i);
+       fHistTrEmcPhiEta[i] = new TH2F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02);
+       fHistTrEmcPhiEta[i]->GetXaxis()->SetTitle("#eta");
+       fHistTrEmcPhiEta[i]->GetYaxis()->SetTitle("#phi");
+       fOutput->Add(fHistTrEmcPhiEta[i]);
+       
+       histname = Form("fHistTrEmcPt_%d",i);
+       fHistTrEmcPt[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
+       fHistTrEmcPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       fHistTrEmcPt[i]->GetYaxis()->SetTitle("counts");
+       fOutput->Add(fHistTrEmcPt[i]);
+       
+       histname = Form("fHistTrPhiEtaNonProp_%d",i);
+       fHistTrPhiEtaNonProp[i] = new TH2F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02);
+       fHistTrPhiEtaNonProp[i]->GetXaxis()->SetTitle("#eta");
+       fHistTrPhiEtaNonProp[i]->GetYaxis()->SetTitle("#phi");
+       fOutput->Add(fHistTrPhiEtaNonProp[i]);
+
+       histname = Form("fHistTrPtNonProp_%d",i);
+       fHistTrPtNonProp[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
+       fHistTrPtNonProp[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       fHistTrPtNonProp[i]->GetYaxis()->SetTitle("counts");
+       fOutput->Add(fHistTrPtNonProp[i]);
+       
+       histname = Form("fHistDeltaEtaPt_%d",i);
+       fHistDeltaEtaPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
+       fHistDeltaEtaPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       fHistDeltaEtaPt[i]->GetYaxis()->SetTitle("#delta#eta");
+       fOutput->Add(fHistDeltaEtaPt[i]);
+       
+       histname = Form("fHistDeltaPhiPt_%d",i);
+       fHistDeltaPhiPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, 200, -2, 2);
+       fHistDeltaPhiPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       fHistDeltaPhiPt[i]->GetYaxis()->SetTitle("#delta#phi");
+       fOutput->Add(fHistDeltaPhiPt[i]);
+       
+       histname = Form("fHistDeltaPtvsPt_%d",i);
+       fHistDeltaPtvsPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, fNbins, -fMaxBinPt/2, fMaxBinPt/2);
+       fHistDeltaPtvsPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       fHistDeltaPtvsPt[i]->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+       fHistDeltaPtvsPt[i]->GetZaxis()->SetTitle("counts");
+       fOutput->Add(fHistDeltaPtvsPt[i]);
+      }
+    }
   }
 
-  if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-    fHistClusPhiEtaEnergy = new TH3F("fHistClusPhiEtaEnergy","Phi-Eta-Energy distribution of clusters", fNbins, fMinBinPt, fMaxBinPt, 80, -2, 2, 128, 0, 6.4);
-    fHistClusPhiEtaEnergy->GetXaxis()->SetTitle("E [GeV]");
-    fHistClusPhiEtaEnergy->GetYaxis()->SetTitle("#eta");
-    fHistClusPhiEtaEnergy->GetZaxis()->SetTitle("#phi");
-    fOutput->Add(fHistClusPhiEtaEnergy);
+  if (fClusterCollArray.GetEntriesFast()>0) {
+    TString histname;
+
+    for (Int_t i = 0; i < fNcentBins; i++) {
+      histname = "fHistClusPhiEtaEnergy_";
+      histname += i;
+      fHistClusPhiEtaEnergy[i] = new TH3F(histname, histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02, fNbins, fMinBinPt, fMaxBinPt);
+      fHistClusPhiEtaEnergy[i]->GetXaxis()->SetTitle("#eta");
+      fHistClusPhiEtaEnergy[i]->GetYaxis()->SetTitle("#phi");
+      fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("E_{cluster} (GeV)");
+      fOutput->Add(fHistClusPhiEtaEnergy[i]);
+
+      histname = "fHistClusDeltaPhiEPEnergy_";
+      histname += i;
+      fHistClusDeltaPhiEPEnergy[i] = new TH2F(histname, histname, fNbins, fMinBinPt, fMaxBinPt, 100, 0, TMath::Pi());
+      fHistClusDeltaPhiEPEnergy[i]->GetXaxis()->SetTitle("E_{cluster} (GeV)");
+      fHistClusDeltaPhiEPEnergy[i]->GetYaxis()->SetTitle("#phi_{cluster} - #psi_{RP}");
+      fOutput->Add(fHistClusDeltaPhiEPEnergy[i]);
+
+      if (fIsEmbedded) {
+       histname = "fHistClusMCEnergyFraction_";
+       histname += i;
+       fHistClusMCEnergyFraction[i] = new TH1F(histname, histname, fNbins, 0, 1.2);
+       fHistClusMCEnergyFraction[i]->GetXaxis()->SetTitle("MC fraction");
+       fHistClusMCEnergyFraction[i]->GetYaxis()->SetTitle("counts");
+       fOutput->Add(fHistClusMCEnergyFraction[i]);
+      }
+    }
+
+    fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins,  -1e-6, 1e-6);
+    fHistClusTimeEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
+    fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
+    fOutput->Add(fHistClusTimeEnergy);
 
-    fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
-    fHistNCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
+    Int_t nbins = fMaxCellsInCluster;
+    while (nbins > fNbins) nbins /= 2;
+    fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, nbins, -0.5, fMaxCellsInCluster-0.5);
+    fHistNCellsEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
     fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
     fOutput->Add(fHistNCellsEnergy);
-  }
 
-  if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-   
-    fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinBinPt, fMaxBinPt);
-    fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
-    fHistCellsEnergy->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistCellsEnergy);
+    fHistFcrossEnergy = new TH2F("fHistFcrossEnergy","fHistFcrossEnergy", fNbins, fMinBinPt, fMaxBinPt, 200, -3.5, 1.5);
+    fHistFcrossEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
+    fHistFcrossEnergy->GetYaxis()->SetTitle("F_{cross}");
+    fOutput->Add(fHistFcrossEnergy); 
+
+    fHistClusEnergyMinusCellEnergy = new TH2F("fHistClusEnergyMinusCellEnergy","fHistClusEnergyMinusCellEnergy", 
+                                             fNbins, fMinBinPt, fMaxBinPt, fNbins, -fMaxBinPt/2, fMaxBinPt/2);
+    fHistClusEnergyMinusCellEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
+    fHistClusEnergyMinusCellEnergy->GetYaxis()->SetTitle("E_{cluster} - #Sigma_{i}E_{cell,i} (GeV)");
+    fOutput->Add(fHistClusEnergyMinusCellEnergy); 
+     
+    fHistCellsAbsIdEnergy = new TH2F("fHistCellsAbsIdEnergy","fHistCellsAbsIdEnergy", 11600,0,11599,(Int_t)(fNbins / 2), fMinBinPt, fMaxBinPt / 2);
+    fHistCellsAbsIdEnergy->GetXaxis()->SetTitle("cell abs. Id");
+    fHistCellsAbsIdEnergy->GetYaxis()->SetTitle("E_{cluster} (GeV)");
+    fHistCellsAbsIdEnergy->GetZaxis()->SetTitle("counts");    
+    fOutput->Add(fHistCellsAbsIdEnergy);
     
     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(), 128, 0, 6.4);
-    fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
-    fOutput->Add(fHistTrackPhi[i]);
-
-    TString histnameeta("fHistTrackEta_");
-    histnameeta += i;
-    fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -2, 2);
-    fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
-    fOutput->Add(fHistTrackEta[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()) {
+       
+  if (fJetCollArray.GetEntriesFast()>0) {
 
     TString histname;
 
-    for (Int_t i = 0; i < 4; i++) {
+    for (Int_t i = 0; i < fNcentBins; i++) {
       histname = "fHistJetsPhiEta_";
       histname += i;
-      fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 80, -2, 2, 128, 0, 6.4);
+      fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1, 1, 101, 0, TMath::Pi() * 2.02);
       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
-      fHistJetsPhiEta[i]->GetZaxis()->SetTitle("p_{T} [GeV/c]");
+      fHistJetsPhiEta[i]->GetZaxis()->SetTitle("counts");
       fOutput->Add(fHistJetsPhiEta[i]);
 
-      histname = "fHistJetsPtNonBias_";
+      histname = "fHistJetsPtArea_";
       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]);
+      fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 50, 0, 1.5);
+      fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+      fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
+      fOutput->Add(fHistJetsPtArea[i]);
+    }
+  }
+  
+  Int_t dim = 0;
+  TString title[20];
+  Int_t nbins[20] = {0};
+  Double_t min[20] = {0};
+  Double_t max[20] = {0};
+  
+  if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
+    title[dim] = "Centrality %";
+    nbins[dim] = 101;
+    min[dim] = 0;
+    max[dim] = 101;
+    dim++;
+
+    if (!fCentMethod2.IsNull()) {
+      title[dim] = Form("Centrality %s %%", fCentMethod2.Data());
+      nbins[dim] = 101;
+      min[dim] = 0;
+      max[dim] = 101;
+      dim++;
+    }
 
-      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 (!fCentMethod3.IsNull()) {
+      title[dim] = Form("Centrality %s %%", fCentMethod3.Data());
+      nbins[dim] = 101;
+      min[dim] = 0;
+      max[dim] = 101;
+      dim++;
+    }
+    
+    if (fDoV0QA==1) {
+      title[dim] = "V0A total multiplicity";
+      nbins[dim] = 200;
+      min[dim] = 0;
+      max[dim] = 20000;
+      dim++;
+
+      title[dim] = "V0C total multiplicity";
+      nbins[dim] = 200;
+      min[dim] = 0;
+      max[dim] = 20000;
+      dim++;
+    }
+    else if (fDoV0QA==2) {
+      title[dim] = "V0A+V0C total multiplicity";
+      nbins[dim] = 300;
+      min[dim] = 0;
+      max[dim] = 30000;
+      dim++;
+    }
 
-      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]);
-      }
+    if (!fRhoName.IsNull()) {
+      title[dim] = "#rho (GeV/c)";
+      nbins[dim] = fNbins*4;
+      min[dim] = 0;
+      max[dim] = 400;
+      dim++;
+    }
 
-      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]);
+    if (fDoEPQA) {
+      title[dim] = "#psi_{RP}";
+      nbins[dim] = 200;
+      min[dim] = -TMath::Pi();
+      max[dim] = TMath::Pi();
+      dim++;
+    }
+  }
 
-      histname = "fHistJetsPtAreaNonBias_";
-      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]);
+  if (fParticleCollArray.GetEntriesFast()>0) {
+    title[dim] = "No. of tracks";
+    nbins[dim] = 3000;
+    min[dim] = -0.5;
+    max[dim] = 6000-0.5;
+    dim++;
+
+    title[dim] = "p_{T,track}^{leading} (GeV/c)";
+    nbins[dim] = fNbins;
+    min[dim] = fMinBinPt;
+    max[dim] = fMaxBinPt;
+    dim++;
+
+    if (fDoLeadingObjectPosition) {
+      title[dim] = "#eta_{track}^{leading}";
+      nbins[dim] = 100;
+      min[dim] = -1;
+      max[dim] = 1;
+      dim++;
+
+      title[dim] = "#phi_{track}^{leading}";
+      nbins[dim] = 101;
+      min[dim] = 0;
+      max[dim] = TMath::Pi() * 2.02;
+      dim++;
+    }
+  }
 
-      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]->GetYaxis()->SetTitle("area");
-      fOutput->Add(fHistJetsPtArea[i]);
+  if (fClusterCollArray.GetEntriesFast()>0) {
+    title[dim] = "No. of clusters";
+    nbins[dim] = 2000;
+    min[dim] = 0;
+    max[dim] = 4000-0.5;
+    dim++;
+
+    title[dim] = "E_{cluster}^{leading} (GeV)";
+    nbins[dim] = fNbins;
+    min[dim] = fMinBinPt;
+    max[dim] = fMaxBinPt;
+    dim++;
+
+    if (fDoLeadingObjectPosition) {
+      title[dim] = "#eta_{cluster}^{leading}";
+      nbins[dim] = 100;
+      min[dim] = -1;
+      max[dim] = 1;
+      dim++;
+
+      title[dim] = "#phi_{cluster}^{leading}";
+      nbins[dim] = 101;
+      min[dim] = 0;
+      max[dim] = TMath::Pi() * 2.02;
+      dim++;
     }
   }
 
+  if (!fCaloCellsName.IsNull()) {
+    title[dim] = "No. of cells";
+    nbins[dim] = 3000;
+    min[dim] = 0;
+    max[dim] = 6000-0.5;
+    dim++;
+  }
+
+  if (fJetCollArray.GetEntriesFast()>0) {
+    title[dim] = "No. of jets";
+    nbins[dim] = 200;
+    min[dim] = 0;
+    max[dim] = 200-0.5;
+    dim++;
+
+    title[dim] = "p_{T,jet}^{leading} (GeV/c)";
+    nbins[dim] = fNbins;
+    min[dim] = fMinBinPt;
+    max[dim] = fMaxBinPt;
+    dim++;
+
+    if (fDoLeadingObjectPosition) {
+      title[dim] = "#eta_{jet}^{leading}";
+      nbins[dim] = 100;
+      min[dim] = -1;
+      max[dim] = 1;
+      dim++;
+
+      title[dim] = "#phi_{jet}^{leading}";
+      nbins[dim] = 101;
+      min[dim] = 0;
+      max[dim] = TMath::Pi() * 2.02;
+      dim++;
+    }
+  }
+
+  fHistEventQA = new THnSparseF("fHistEventQA","fHistEventQA",dim,nbins,min,max);
+  for (Int_t i = 0; i < dim; i++)
+    fHistEventQA->GetAxis(i)->SetTitle(title[i]);
+  fOutput->Add(fHistEventQA);
+
   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 }
 
+//________________________________________________________________________
+void AliAnalysisTaskSAQA::ExecOnce()
+{
+  AliAnalysisTaskEmcalJet::ExecOnce();
+  
+  if (fDoV0QA) {
+    fVZERO = InputEvent()->GetVZEROData();
+    if (!fVZERO) {
+      AliError("AliVVZERO not available");
+    }
+  }
+}
+
 //________________________________________________________________________
 Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
 {
@@ -382,82 +528,174 @@ Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
     return kFALSE;
 
-  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;
+  if (!fCentMethod2.IsNull() || !fCentMethod3.IsNull()) {
+    if (fBeamType == kAA || fBeamType == kpA ) {
+      AliCentrality *aliCent = InputEvent()->GetCentrality();
+      if (aliCent) {
+       if (!fCentMethod2.IsNull()) 
+         fCent2 = aliCent->GetCentralityPercentile(fCentMethod2); 
+       if (!fCentMethod3.IsNull()) 
+         fCent3 = aliCent->GetCentralityPercentile(fCentMethod3);
       }
     }
   }
 
+  if (fVZERO) {
+    fV0ATotMult = AliESDUtils::GetCorrV0A(fVZERO->GetMTotV0A(),fVertex[2]);
+    fV0CTotMult = AliESDUtils::GetCorrV0C(fVZERO->GetMTotV0C(),fVertex[2]);
+  }
+
   return kTRUE;
 }
 
+
 //________________________________________________________________________
 Bool_t AliAnalysisTaskSAQA::FillHistograms()
 {
   // Fill histograms.
 
-  fHistCentrality->Fill(fCent);
-  if (fTracks)
-    fHistTracksCent->Fill(fCent, fTracks->GetEntriesFast());
-  if (fCaloClusters)
-    fHistClusCent->Fill(fCent, fCaloClusters->GetEntriesFast());
+  Float_t trackSum = 0;
+  Float_t clusSum = 0;
+  Float_t cellSum = 0;
+  Float_t cellCutSum = 0;
 
-  fHistZVertex->Fill(fVertex[2]);
+  Int_t ntracks = 0;
+  Int_t nclusters = 0;
+  Int_t ncells = 0;
+  Int_t njets = 0;
 
-  Float_t trackSum = DoTrackLoop();
+  Float_t leadingClusE = 0;
+  Float_t leadingClusEta = 0;
+  Float_t leadingClusPhi = 0;
 
-  DoJetLoop();
+  Float_t leadingTrackPt = 0;
+  Float_t leadingTrackEta = 0;
+  Float_t leadingTrackPhi = 0;
 
-  if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
+  Float_t leadingJetPt = 0;
+  Float_t leadingJetEta = 0;
+  Float_t leadingJetPhi = 0;
+    
+  if (fTracks) {
+    AliVParticle *leadingTrack = 0;
 
-    if (fTracks && fCaloClusters)
-      fHistClusTracks->Fill(fTracks->GetEntriesFast(), fCaloClusters->GetEntriesFast());
+    ntracks = DoTrackLoop(trackSum, leadingTrack);
+    AliDebug(2,Form("%d tracks found in the event", ntracks));
 
-    Float_t clusSum = DoClusterLoop();
+    if (leadingTrack) {
+      leadingTrackPt = leadingTrack->Pt();
+      leadingTrackEta = leadingTrack->Eta();
+      leadingTrackPhi = leadingTrack->Phi();
+    }
+  } 
 
-    Float_t cellSum = 0, cellCutSum = 0;
-    
-    Int_t ncells = DoCellLoop(cellSum, cellCutSum);
+  if (fCaloClusters) {
+    AliVCluster  *leadingClus = 0;
+
+    nclusters = DoClusterLoop(clusSum, leadingClus);
+    AliDebug(2,Form("%d clusters found in the event", nclusters));
 
-    if (fTracks)
-      fHistCellsTracks->Fill(fTracks->GetEntriesFast(), ncells);
+    if (leadingClus) {
+      TLorentzVector leadingClusVect;
+      leadingClus->GetMomentum(leadingClusVect, fVertex);
+      leadingClusE = leadingClus->E();
+      leadingClusEta = leadingClusVect.Eta();
+      leadingClusPhi = leadingClusVect.Phi();
+    }
 
-    fHistCellsCent->Fill(fCent, ncells);
+    fHistChVSneClus->Fill(clusSum, trackSum);
+  }
+  
+  if (fCaloCells) {
+    ncells = DoCellLoop(cellSum, cellCutSum);
+    AliDebug(2,Form("%d cells found in the event", ncells));
     
     fHistChVSneCells->Fill(cellSum, trackSum);
-    fHistChVSneClus->Fill(clusSum, trackSum);
     fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
+  }
 
-    if (fDoTrigger) {
-      Float_t maxTrgClus = DoTriggerClusLoop();
-      fHistMaxL1ClusCent->Fill(fCent, maxTrgClus);
-    
-      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);
+  if (fJets) {
+    AliEmcalJet  *leadingJet = 0;
+
+    njets = DoJetLoop(leadingJet);
+    AliDebug(2,Form("%d jets found in the event", njets));
+
+    if (leadingJet) {
+      leadingJetPt = leadingJet->Pt();
+      leadingJetEta = leadingJet->Eta();
+      leadingJetPhi = leadingJet->Phi();
     }
   }
 
+  FillEventQAHisto(fCent, fCent2, fCent3, fV0ATotMult, fV0CTotMult, fEPV0, fRhoVal, 
+                  ntracks, nclusters, ncells, njets, 
+                  leadingTrackPt, leadingTrackEta, leadingTrackPhi, 
+                  leadingClusE, leadingClusEta, leadingClusPhi, 
+                  leadingJetPt, leadingJetEta, leadingJetPhi);
+
   return kTRUE;
 }
 
+//________________________________________________________________________
+void AliAnalysisTaskSAQA::FillEventQAHisto(Float_t cent, Float_t cent2, Float_t cent3, Float_t v0a, Float_t v0c, 
+                                          Float_t ep, Float_t rho, Int_t ntracks, Int_t nclusters, Int_t ncells, Int_t njets, 
+                                          Float_t maxTrackPt, Float_t maxTrackEta, Float_t maxTrackPhi,
+                                          Float_t maxClusterE, Float_t maxClusterEta, Float_t maxClusterPhi,
+                                          Float_t maxJetPt, Float_t maxJetEta, Float_t maxJetPhi)
+{
+  Double_t contents[20]={0};
+
+  for (Int_t i = 0; i < fHistEventQA->GetNdimensions(); i++) {
+    TString title(fHistEventQA->GetAxis(i)->GetTitle());
+    if (title=="Centrality %")
+      contents[i] = cent;
+    else if (title==Form("Centrality %s %%", fCentMethod2.Data()))
+      contents[i] = cent2;
+    else if (title==Form("Centrality %s %%", fCentMethod3.Data()))
+      contents[i] = cent3;
+    else if (title=="V0A total multiplicity")
+      contents[i] = v0a;
+    else if (title=="V0C total multiplicity")
+      contents[i] = v0c;
+    else if (title=="V0A+V0C total multiplicity")
+      contents[i] = v0a+v0c;
+    else if (title=="#psi_{RP}")
+      contents[i] = ep;
+    else if (title=="#rho (GeV/c)")
+      contents[i] = rho;
+    else if (title=="No. of tracks")
+       contents[i] = ntracks;
+    else if (title=="No. of clusters")
+      contents[i] = nclusters;
+    else if (title=="No. of cells")
+      contents[i] = ncells;
+    else if (title=="No. of jets")
+      contents[i] = njets;
+    else if (title=="p_{T,track}^{leading} (GeV/c)")
+      contents[i] = maxTrackPt;
+    else if (title=="#eta_{track}^{leading}")
+      contents[i] = maxTrackEta;
+    else if (title=="#phi_{track}^{leading}")
+      contents[i] = maxTrackPhi;
+    else if (title=="E_{cluster}^{leading} (GeV)")
+      contents[i] = maxClusterE;
+    else if (title=="#eta_{cluster}^{leading}")
+      contents[i] = maxClusterEta;
+    else if (title=="#phi_{cluster}^{leading}")
+      contents[i] = maxClusterPhi;
+    else if (title=="p_{T,jet}^{leading} (GeV/c)")
+      contents[i] = maxJetPt;
+    else if (title=="#eta_{jet}^{leading}")
+      contents[i] = maxJetEta;
+    else if (title=="#phi_{jet}^{leading}")
+      contents[i] = maxJetPhi;
+    else 
+      AliWarning(Form("Unable to fill dimension %s!",title.Data()));
+  }
+
+  fHistEventQA->Fill(contents);
+}
+
 //________________________________________________________________________
 Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
 {
@@ -471,8 +709,9 @@ Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
   const Int_t ncells = cells->GetNumberOfCells();
 
   for (Int_t pos = 0; pos < ncells; pos++) {
-    Float_t amp = cells->GetAmplitude(pos);
-    fHistCellsEnergy->Fill(amp);
+    Float_t amp   = cells->GetAmplitude(pos);
+    Int_t   absId = cells->GetCellNumber(pos);
+    fHistCellsAbsIdEnergy->Fill(absId,amp);
     sum += amp;
     if (amp < fCellEnergyCut)
       continue;
@@ -483,17 +722,90 @@ Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
 }
 
 //________________________________________________________________________
-Float_t AliAnalysisTaskSAQA::DoClusterLoop()
+Double_t AliAnalysisTaskSAQA::GetCellEnergySum(AliVCluster *cluster, AliVCaloCells *cells)
+{
+  Double_t sum = 0;
+  for (Int_t i = 0; i < cluster->GetNCells(); i++) 
+    sum += cells->GetCellAmplitude(cluster->GetCellAbsId(i));
+  return sum;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskSAQA::GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
+{
+  Int_t    AbsIdseed  = -1;
+  Double_t Eseed      = 0;
+  for (Int_t i = 0; i < cluster->GetNCells(); i++) {
+    if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > AbsIdseed) {
+      Eseed     = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
+      AbsIdseed = cluster->GetCellAbsId(i);
+    }
+  }
+
+  if (Eseed < 1e-9)
+    return 100;
+
+  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
+  fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta); 
+  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);  
+  
+  //Get close cells index and energy, not in corners
+  
+  Int_t absID1 = -1;
+  Int_t absID2 = -1;
+  
+  if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
+  if (iphi > 0)                                 absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
+  
+  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
+  
+  Int_t absID3 = -1;
+  Int_t absID4 = -1;
+  
+  if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
+    absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
+    absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta-1); 
+  }
+  else if (ieta == 0 && imod%2) {
+    absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta+1);
+    absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1); 
+  }
+  else {
+    if (ieta < AliEMCALGeoParams::fgkEMCALCols-1) 
+      absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
+    if (ieta > 0)                                 
+      absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1); 
+  }
+  
+  Double_t  ecell1 = cells->GetCellAmplitude(absID1);
+  Double_t  ecell2 = cells->GetCellAmplitude(absID2);
+  Double_t  ecell3 = cells->GetCellAmplitude(absID3);
+  Double_t  ecell4 = cells->GetCellAmplitude(absID4);
+
+  Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
+  
+  Double_t Fcross = 1 - Ecross/Eseed;
+
+  return Fcross;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskSAQA::DoClusterLoop(Float_t &sum, AliVCluster* &leading)
 {
   // Do cluster loop.
 
   if (!fCaloClusters)
     return 0;
 
-  Float_t sum = 0;
+  Int_t nAccClusters = 0;
+
+  AliVCaloCells *cells = InputEvent()->GetEMCALCells();
+
+  sum = 0;
+  leading = 0;
 
   // Cluster loop
-  Int_t nclusters =  fCaloClusters->GetEntriesFast();
+  const Int_t nclusters = fCaloClusters->GetEntriesFast();
 
   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
@@ -502,37 +814,56 @@ Float_t AliAnalysisTaskSAQA::DoClusterLoop()
       continue;
     }  
 
-    if (!AcceptCluster(cluster, kTRUE))
+    if (!AcceptCluster(cluster))
       continue;
 
     sum += cluster->E();
 
+    if (!leading || leading->E() < cluster->E()) leading = cluster;
+
     TLorentzVector nPart;
     cluster->GetMomentum(nPart, fVertex);
 
-    fHistClusPhiEtaEnergy->Fill(cluster->E(), nPart.Eta(), nPart.Phi());
+    fHistClusPhiEtaEnergy[fCentBin]->Fill(nPart.Eta(), nPart.Phi(), cluster->E());
+
+    Double_t ep = nPart.Phi() - fEPV0;
+    while (ep < 0) ep += TMath::Pi();
+    while (ep >= TMath::Pi()) ep -= TMath::Pi();
+    fHistClusDeltaPhiEPEnergy[fCentBin]->Fill(cluster->E(), ep);
+
     fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
 
     fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
+
+    if (cells) fHistFcrossEnergy->Fill(cluster->E(), GetFcross(cluster, cells));
+
+    if (cells) fHistClusEnergyMinusCellEnergy->Fill(cluster->E(), cluster->E() - GetCellEnergySum(cluster,cells));
+
+    if (fHistClusMCEnergyFraction[fCentBin])
+      fHistClusMCEnergyFraction[fCentBin]->Fill(cluster->GetMCEnergyFraction());
+
+    nAccClusters++;
   }
 
-  return sum;
+  return nAccClusters;
 }
 
 //________________________________________________________________________
-Float_t AliAnalysisTaskSAQA::DoTrackLoop()
+Int_t AliAnalysisTaskSAQA::DoTrackLoop(Float_t &sum, AliVParticle* &leading)
 {
   // Do track loop.
 
   if (!fTracks)
     return 0;
 
-  Float_t sum = 0;
+  Int_t nAccTracks = 0;
 
-  Int_t ntracks = fTracks->GetEntriesFast();
-  Int_t nclusters = 0;
-  if (fCaloClusters)
-    nclusters = fCaloClusters->GetEntriesFast();
+  sum = 0;
+  leading = 0;
+
+  const Int_t ntracks = fTracks->GetEntriesFast();
+  Int_t neg = 0;
+  Int_t zero = 0;
 
   for (Int_t i = 0; i < ntracks; i++) {
 
@@ -543,97 +874,85 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
       continue; 
     }
 
-    AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
-    
-    if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
+    if (!AcceptTrack(track)) 
       continue;
-    
-    fHistTracksPt->Fill(track->Pt());
+
+    nAccTracks++;
 
     sum += track->P();
 
-    Int_t label = track->GetLabel();
-      
-    fHistTrPhiEta->Fill(track->Eta(), track->Phi());
-    
-    fHistTrackEta[4]->Fill(track->Eta());
-    fHistTrackPhi[4]->Fill(track->Phi());
+    if (!leading || leading->Pt() < track->Pt()) leading = track;
 
-    if (label >= 0 && label < 4) {
-      fHistTrackEta[label]->Fill(track->Eta());
-      fHistTrackPhi[label]->Fill(track->Phi());
+    if (fParticleLevel) {
+      fHistTrPhiEtaPt[fCentBin][0]->Fill(track->Eta(), track->Phi(), track->Pt());
     }
+    else {
+      fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
+      if (track->GetLabel() == 0) {
+       zero++;
+       if (fHistTrPhiEtaZeroLab[fCentBin]) {
+         fHistTrPhiEtaZeroLab[fCentBin]->Fill(track->Eta(), track->Phi());
+         fHistTrPtZeroLab[fCentBin]->Fill(track->Pt());
+       }
+      }
 
-    if (!vtrack)
-      continue;
+      if (track->GetLabel() < 0)
+       neg++;
 
-    if (vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999)
-      fHistTrPhiEtaNonProp->Fill(vtrack->Eta(), vtrack->Phi());
+      Int_t type = 0;
 
-    fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
-    fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
-    fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
+      AliPicoTrack* ptrack = dynamic_cast<AliPicoTrack*>(track);
+      if (ptrack)
+       type = ptrack->GetTrackType();
 
-    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());
+      if (type >= 0 && type < 3)
+       fHistTrPhiEtaPt[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
+      else
+       AliDebug(2,Form("%s: track type %d not recognized!", GetName(), type));
     }
-  }
-  
-  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);
+    AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
+
+    if (!vtrack)
+      continue;
+
+    if ((vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999) && fHistTrPhiEtaNonProp[fCentBin]) {
+      fHistTrPhiEtaNonProp[fCentBin]->Fill(vtrack->Eta(), vtrack->Phi());
+      fHistTrPtNonProp[fCentBin]->Fill(vtrack->Pt());
     }
+
+    if (fHistTrEmcPhiEta[fCentBin])
+      fHistTrEmcPhiEta[fCentBin]->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());   
+    if (fHistTrEmcPt[fCentBin])
+      fHistTrEmcPt[fCentBin]->Fill(vtrack->GetTrackPtOnEMCal());   
+    if (fHistDeltaEtaPt[fCentBin])
+      fHistDeltaEtaPt[fCentBin]->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
+    if (fHistDeltaPhiPt[fCentBin])
+      fHistDeltaPhiPt[fCentBin]->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
+    if (fHistDeltaPtvsPt[fCentBin])
+      fHistDeltaPtvsPt[fCentBin]->Fill(vtrack->Pt(), vtrack->Pt() - vtrack->GetTrackPtOnEMCal());
   }
-    
-  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();
+
+  if (fHistTrNegativeLabels)
+    fHistTrNegativeLabels->Fill(1. * neg / ntracks);
+
+  if (fHistTrZeroLabels)
+    fHistTrZeroLabels->Fill(1. * zero / ntracks);
+
+  return nAccTracks;
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskSAQA::DoJetLoop()
+Int_t AliAnalysisTaskSAQA::DoJetLoop(AliEmcalJet* &leading)
 {
   // Do jet loop.
 
   if (!fJets)
-    return;
+    return 0;
+
+  Int_t nAccJets = 0;
+
+  leading = 0;
 
   Int_t njets = fJets->GetEntriesFast();
 
@@ -646,84 +965,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());
-
-    if (jet->MaxTrackPt() > fPtBiasJetTrack)
-      fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
-    
-    if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus)
-      fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
-    
-    if (!AcceptBiasJet(jet))
-      continue;
+    if (!leading || leading->Pt() < jet->Pt()) leading = jet;
 
-    fHistJetsPt[fCentBin]->Fill(jet->Pt());
-    fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
+    nAccJets++;
 
     fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
+    fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
   }
-}
-
-//________________________________________________________________________
-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 *) 
-{
-  // Called once at the end of the analysis.
+  return nAccJets;
 }